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ABSTRACT 


A  numerical  study  is  conducted  to  study  lithium-sulfur  hexafluoride  (Li-SF6) 
wick  diffusion  flames.  The  objective  of  this  study  is  to  assess  the  effects  of  changing  the 
geometry  (height  H  and  aspect  ratio  H/yo)  and  ambient  conditions  (free-stream  velocity 
Uoo  and  gravity)  on  the  burning  rate  and  heat  transfer.  Wick  combustion  is  identified  as  a 
boundary-layer  gaseous  diffusion  flame  with  multiphase  combustion  products.  A 
mathematical  model  for  wick  diffusion  flames  is  established  employing  a  conserved 
scalar  approach.  Both  forced  and  mixed  convective  burning  conditions  are  considered. 
Laminar,  variable-property,  boundary-layer  equations  are  cast  into  dimensionless  forms 
using  the  modified  Howarth-Dorodnitzyn  transformation  for  vertical  plates  and  cylinders. 
The  state  relationships  for  the  properties  are  taken  from  existing  data  of  Li-SF6 
combustion  at  a  pressure  of  P  =  0.01  MPa. 

Forced  convective  burning,  for  which  the  Reynolds  number  (Re  =  UeoHAoc)  is  the 
important  dimensionless  parameter,  is  studied  first.  The  results  show  that  increasing  Uoo 
increases  the  total  burning  rate,  rh.  A  relationship  between  m  and  Re  is  obtained  for  both 
planar  and  cylindrical  wicks  of  a  given  geometry  (H  =  100  mm,  yo  =  12.5  mm).  The  flat- 
plate  solution  yields  a  fuel  mass  burning  rate  per  unit  surface  area  (i.e.  local  fuel  burning 
rate)  following  the  x’^/2  dependence  of  the  classical  similarity  solution,  where  x  is  the 
streamwise  distance.  This  dependence  is  not  evident  in  the  cylinder  solution.  Cylindrical 
wick  geometries  yield  enhanced  burning  rates  over  planar  wicks. 

For  the  case  of  mixed  convective  burning,  the  burning  rate  results  approach  either 
the  forced  or  natural  convective  burning  limits  as  u«  is  increased  or  decreased, 
respectively.  Critical  Richardson  numbers  (Ri  =  Gr/Re^)  specifying  these  burning  limits 
are  determined  for  the  given  baseline  condition.  When  the  wick  height  is  the  only 
parameter  varied,  the  different  curves  for  the  local  fuel  burning  rate,m",  collapse  on  one 
another,  if  plotted  versus  the  local  Richardson  number  (Rix).  Reducing  gravity  results  in 
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a  lower  burning  rate  because  the  influence  of  natural  convection  is  diminished.  Under 
microgravity  conditions  (1/1 000th  of  gravity  at  sea  level),  mixed  convective  burning 
nearly  resembles  forced  convection.  The  results  for  cylindrical  wicks  approach  the 
results  for  planar  wicks  as  the  radius  is  increased,  keeping  the  height  constant.  The  local 
fuel  burning  rates  curves  for  cylindrical  wicks  of  50.8  and  100  mm  (H/yo  =  2  and  1, 
respectively)  deviate  slightly  from  the  curve  for  a  flat  plate  undergoing  the  same  burning 
conditions.  In  general,  cylindrical  wick  geometries  yield  enhanced  burning  rates  over 
planar  wicks. 

The  major  limitation  of  the  present  study  is  that  the  analysis  does  not  take  in 
account  radiative  heat  transfer,  end  effects  of  a  finite  wick,  and  the  presence  of  non¬ 
condensable  gases.  Further  analysis  to  account  for  these  effects  is  recommended. 
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CHAPTER  1 
INTRODUCTION 


A  novel  liquid-metal  combustion  system  involving  the  reaction  of  lithium  (Li)  and 
sulfur  hexafluoride  (SF6)  is  numerically  studied.  This  system  (e.g.,  see  Groff  and  Faeth, 
1978)  is  unique  in  that  it  has  a  high  energy  density  and  a  high  specific  energy  and  that  it 
yields  condensed-phase  products  under  normal  operating  conditions.  The  high  energy 
density  and  high  specific  energy  satisfy  the  volume  and  weight  requirements  of 
propulsion  applications  and  the  condensed-phase  product  provides  a  means  of  closed 
system  operation  as  desired  by  deep-sea  propulsion  vessels;  thus  the  liquid-metal 
combustion  system  is  ideal  for  undersea  applications. 

The  present  study  is  a  one- year  continuation  of  the  research  of  Chen  et  al.  (1991). 
This  report  summarizes  the  research  accomplished  following  that  repon;  namely,  on  the 
computational  study  of  forced  convective  burning  of  Li-SFs  wick  diffusion  flames. 

This  Chapter  summarizes  the  literature  review  and  the  objectives  of  the  study. 
The  publications  resulting  from  this  Navy  sponsored  research,  i.e.,  Chen  et  al.  (1991)  and 
the  present  study,  include  Chen  et  al.  (1990  and  1991),  Hsu  and  Chen  (1991  and  1992), 
Lyu  and  Chen  (1991,  1991a  and  1991b),  Wu  and  Chen  (1992),  Damaso  and  Chen  (1992). 
Two  doctoral  dissertations,  Lyu  (1991)  and  Hsu  (1991),  were  supported  by  this  Navy 
sponsored  research  program;  and  one  additional  doctoral  dissertation,  Wu  (1991),  and 
one  M.  S.  thesis,  Damaso  (1992),  were  partially  supported. 

Li-SF^  Combustion 


The  reaction  of  Li  and  SFg  follows  the  stoichiometry: 
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8  Li  +  SF5  — >  6  LiF  +  Li2S  ( 1  - 1 ) 

The  combustion  products  (LiF  and  Li2S)  are  soluble  in  liquid  lithium,  yielding  two 
immiscible  liquids-  fuel-rich  and  product-rich  liquids.  The  density  of  the  product-rich 
liquid  is  heavier  than  the  fuel-rich  liquid  (and  pure  Li),  allowing  for  constant  volume 
operation.  This  system  is  unique  in  that  it  has  a  high  energy  density  and  a  high  specific 
energy  and  that  it  yields  condensed-phase  products  under  normal  operating  conditions. 
The  high  energy  density  and  high  specific  energy  satisfy  the  volume  and  weight 
requirements  of  undersea  propulsion  systems,  as  well  as  space  applications.  The 
condensed-phase  product  provides  a  means  of  closed  system  operation.  Furthermore,  Li 
is  stable  (non-reacting),  and,  at  room  temperature,  is  in  solid  phase.  SF6  is  non-toxic  and 
stored  as  a  liquid  with  a  high  vapor  pressure.  The  high  vapor  pressure  allows  the 
injection  of  SFg  into  the  combustor  without  the  need  of  a  pressurizing  device.  Li  and  SFe 
thus  make  an  attractive  reaction  pair  for  undersea  propulsion  and  space  applications. 

One  major  system  configuration  already  studied  is  the  reactive  heat  pipe  (or  wick 
configuration),  e.g.  see  Lyu  et  al.  (1990)  and  the  references  cited  therein.  When  a  wick  is 
employed,  liquid  Li  is  supplied  through  capillary  action,  and  the  heat  transfer  is 
accomplished  by  evaporation  and  condensation  of  Li.  While  the  wick  stands  vertically 
during  combustion,  the  buoyancy  effects  result  in  natural  convective  burning  along  the 
surface.  However,  in  space  applications,  where  microgravity  conditions  exist,  buoyancy 
effects  may  be  negligible.  This  leads  to  the  introduction  of  forced  convective  burning,  i” 
which  a  free  stream  of  oxidant  (SF6)  is  blown  past  the  wick  surface  to  facilitate 
combustion.  One  advantage  forced  convective  burning  has  is  that  by  varying  the  free 
stream  velocity  of  the  oxidant  the  fuel  burning  rate  can  be  changed  accordingly,  thus. 
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allowing  for  changes  in  heat  output  of  the  combustor.  Furthermore,  forced  convective 
burning  is  independent  of  gravity. 


Previous  Work 

Previous  work  has  been  done  in  the  area  of  convective  burning  and,  in  panicular, 
natural  convective  burning.  Regarding  natural  convective  burning  of  Li-SF6  along  a  flat 
plate,  Lyu  (1991)  showed  that  the  local  fuel  burning  rate,  m",  is  proportional  to 
where  ^  is  the  dimensionless  streamwise  distance.  In  addition,  Lyu  obtained  similarity 
solutions  for  the  planar  (flat  plate)  wick  geometry,  based  on  a  system  pressure  of  0.01 
MPa  and  wicks  of  height  (H  =  100  mm).  Cylindrical  wick  geometries  with  radius  (yo  = 
12.5  mm)  did  not  yield  similarity  solutions.  Much  of  this  study  is  based  on  Lyu's  work 
and  extends  it  into  the  realm  of  forced  and  mixed  convective  burning.  Regarding  forced 
convective  burning,  Williams  (1985)  predicted  the  local  fuel  burning  rate  along  a 
horizontal  fuel  plate  to  be  proportional  to 

Though  forced  convective  burning  is  ideal  for  space  applications,  experimental 
research  is  presently  confined  to  earth,  where  microgravity  conditions  do  not  exist.  In 
these  experiments,  both  modes  of  burning— natural  and  forced— must  be  accounted  for. 
Therefore,  the  study  of  mixed  convective  burning  gains  its  significance  in  relation  to 
space  applications.  A  major  challenge  is  to  perform  earth-bound  experiments  (involving 
mixed  convective  burning)  that  approximate  conditions  found  in  space.  To  aid  in  these 
experiments,  this  study  presents  results  regarding  mixed  convective  burning. 

This  study  draws  on  the  work  done  for  natural  and  forced  convective  burning,  as 
well  as  for  mixed  convection.  Jaluria  (1986)  did  a  numerical  study  on  aiding  mixed 
convective  flow  (non-reacting)  over  a  heated  vertical  surface.  The  study  concluded  that 
natural  convection  has  a  major  influence  in  mixed  convection.  In  addition,  the  study 
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found  that  for  large  values  of  the  mixed  convection  parameter,  Gr/Re5/2,  the  Nusselt 
number  curves  approached  the  results  for  pure  natural  convection.  Yao  (1985)  studied 
mixed  convection  along  a  vertical  flat  plate.  Yao’s  paper  stated  that  natural  convection 
effects  are  small  near  the  leading  edge,  and  that  the  forced  convection  limit  is  the  solution 
of  the  problem  at  the  leading  edge.  Yao  also  presented  results  for  heat  transfer  in  terms 
of  the  Nusselt  number. 


Objectives 

As  discussed  earlier,  the  wick-type  liquid-metal  combustor  is  of  interest  due  to  its 
use  as  an  energy  source  for  space  applications.  The  overall  objective  of  the  present  study 
is  to  extend  the  methodology  (developed  for  lithium-sulfur  hexafluoride  wick  diffusion 
flames  under  natural  convective  burning  conditions)  to  investigate  forced  convective  and 
mixed  convective  burning  of  wick  diffusion  flames.  Upon  modeling  forced  and  mixed 
convective  burning,  numerical  solutions  are  obtained  to  assess  the  effects  due  to 
geometry  (cylinder  versus  flat  plate,  wick  height,  aspect  ratio)  and  physical  environment 
(free-stream  velocity  of  oxidant,  gravity)  on  the  wick  combustion  of  Li  and  SF6. 
Specifically,  the  objectives  of  the  present  study  are  the  following: 

(1)  To  extend  the  mathematical  model  developed  for  wick-type  liquid-metal  diffusion 
flames  undergoing  natural  convective  burning  to  flames  undergoing  mixed 
convective  burning. 

(2)  To  obtain  numerical  solutions  to  assess  the  effects  due  to  geometry  (cylindrical 
versus  planar  wicks,  wick  height,  aspect  ratio)  and  ambient  conditions  (free- 
stream  velocity  of  oxidant,  gravity)  on  the  wick  combustion  of  Li  and  SFg. 

(3)  To  present  results  that  add  to  the  knowledge  of  earth-bound  combustion 
experiments  performed  for  varying  gravity  conditions. 
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CHAPTER  2 

THEORETICAL  CONSIDERATIONS 

The  present  study  employs  a  conserved  scalar  approach  for  wick  diffusion  flames. 
The  governing  equations  and  boundary  conditions  discussed  in  this  chapter  are  cast  into 
dimensionless  forms  employing  a  similarity  transformation.  The  dimensionless 
parameters  or  similarity  variables  are  helpful  in  identifying  the  physical  quantities 
important  to  the  combustion  process. 
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Figure  2-1.  The  coordinate  system. 


The  following  analysis  extends  the  modeling  of  Lyu  (1991)  to  account  for  purely 
forced  convective  burning  (e.g.  that  at  zero-gravity  condition)  and  mixed  convective 
burning  (combined  forced  and  natural  convection).  Following  Lyu  et  al.  (1990),  Li-SF6 
wick  combustion  is  assumed  as  a  diffusion  flame  with  infinitely  fast  chemical  reactions. 
Steady-state,  two-dimensional  boundary-layer  flows  are  also  assumed.  Introducing  the 
assumptions  of  unity  Lewis  number,  identical  (binary)  diffusion  coefficients  for  all 
species,  no  radiative  heat  less,  negligible  viscous  dissipation,  uniform  interface  condition, 
a  conserved  scalar  approach  can  be  employed.  Such  an  approach,  which  combines  the 
species  and  energy  equations  into  a  single  mixture  fraction  equation,  is  applied  in 
modeling  diffusion  flames  when  locally  kinematic  (velocity),  thermal  (temperature),  and 
phase  equilibria  are  established  in  the  flow.  In  this  study,  mixture  fraction  is  defined  as 
the  effective  mass  fraction  resulting  from  the  fluid  supplied  through  the  wick  action  at 
any  point  in  the  flow  field. 

Considering  the  coordinate  system  illustrated  in  Figure  2-1,  the  governing 
equations  are 

Continuity: 

d(puyn)  a(pvyn) 

ax 

Momentum: 


^  au  .  au  \  d  f  „  au^  ,  . 

(y  P57  ag  80  (p-  -  p) 


(2.2) 
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Mixture  Fraction: 


az  az 


(2.3) 


where  go  is  the  gravitational  acceleration  at  sea  level,  ag  is  an  arbitrary  constant  for 
reduced  gravity  conditions;  e.g.  ag  =  0  for  zero  gravity  (or  pure  forced  convection)  and  ag 
=  1  for  normal  gravity  (or  mixed  convection).  The  superscript  "n"  accounts  for  different 
geometries,  i.e.,  n  =  0  for  flat  plates  and  n  =  1  for  cylinders.  D  is  the  binary  diffusion 
coefficient,  and  Z  is  the  mixture  fraction,  i.e.,  a  Shvab-Zeldovich  variable  (Williams, 
1985).  The  mixture  fraction  can  be  viewed  as  a  dimensionless  enthalpy  (h),  or  a  species 
concentration  (Y); 


Z  = 


h  -  hop  _  Ym  *  Y  moo 
hw  -  hoo  Y  mw  *  Y  moo 


(2.3a) 


Eq.  (2.3)  uses  the  conserved  scalar,  Z,  to  combine  the  species  and  energy  equations  (Lyu, 
1991).  The  equations  for  species  (2.3b)  and  energy  (2.3c)  are 


Species: 


5Ym  3Ym 


Energy: 


(2.3b) 


^  ah  ^  ^  ah  I  d  f  .k  dh\ 


(2.3c) 
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Boundary  conditions  are  specified  for  no-slip  walls.  Still  environments  are 
assigned  for  natural  convection  conditions,  and  parallel  flow  is  assigned  for  forced  and 
mixed  convection. 


o 

II 

X 

y  =  0 

y  — ^  oo 

v  =  Z  =  0. 

u  =  0  (natural  convection) 

u  =  Uoo  (forced  and  mixed  convection) 

Z  =  0 

0<x^H 

O 

II 

II 

N 

> 

II 

> 

O 

II 
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0<x^H 

y  — ^  oo 

u  =  0  (natural  convection) 

u  =  Uoo  (forced  and  mixed  convection) 

Z  =  0  (2.4) 

where  Vw  is  determined  (Lyu,  1991)  by  the  heat  and  mass  balance  at  the  interface,  and  u» 
is  the  free-stream  velocity.  The  formulation  presented,  i.e.,  Eqs.  (2.1)  to  (2.4),  is  for 
laininar,  boundary-layer  reacting  flows. 

Similarity  Transformation 

The  governing  equations  are  cast  into  a  dimensionless  form  in  terms  of  the  stream 
function  and  similarity  variables,  using  a  modified  Howarth-Dorodnitzyn  transformation. 
Tl.e  transformation  has  independent  variables  ^  and  T),  which  are  defined  as 


9 


The  term  d^dy  is  zero  because  4  is  not  a  function  of  y.  Evaluating  the  partial  derivatives, 
one  obtains 
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^  ^  [2^Hpoo|icoUoo]l/2  y" 

3t1  dT|  ° 


(2.10) 


and 


_ ^]l/2o  yH 

L2^poomJ  y" 


(2.11) 


Substituting  Eqs.  (2.9),  (2.10),  and  (2.11)  into  Eq.  (2.6),  an  expression  for  u  is  obtained: 


af 

U  =Uoo  — 

an 


(2.12) 


Similarly,  an  expression  for  v  can  be  obtained  by  first  taking  the  partial  derivative 
of  V  with  respect  to  x: 


av  _  ^  a^  ay  an 


(2.13) 


or 


^  ^  (2^P«MooU«)^^  yj  ^  ^  [f  (2^HpooMcoU«)‘^  y"]  ^ 


(2.14) 


Evaluating  Eq.  (2.14)  and  substituting  into  Eq.  (2.7),  an  expression  for  v  is  obtained: 


1  y"  /^pooPooUoo^/2 


py 


24H 


a^  ^  a^ 


(2.15) 
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With  expressions  for  u  and  v  derived,  Bu/dx  and  Bu/By  must  be  derived  to 
complete  the  evaluation  of  the  left-hand  side  of  the  momentum  equation,  Eq.  (2.3).  First, 
the  term  3u/3x  is  found: 


9u  B  ^  Bf\ 


(2.16) 


9u  2. 
ax"““[H 


a2f  ^  an  92f' 
a^ari  aTi2 


(2.17) 


Then,  the  term  au/ay  is  expressed  as 


au  B  /■  Bf\  a^f  aq 


(2.18) 


^  y"  L2^pooMoJ 


(2.19) 


The  terms  u(au/ax)  and  v(au/ay)  can  now  be  evaluated. 


t7U 


— r 

1  a2f 

aq  J 

^a^aq 

(2.20) 
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3u  Uoo^ 

~ '  2%n 


^  ^^aTiaTi2 


(2.21) 


With  the  terms  on  the  left-hand  side  of  the  momentum  equation  determined,  the 
focus  of  attention  now  turns  to  the  right-hand  side.  Multiplying  Eq.  (2.19)  (the 
expression  for  du/dy)  by  y"  and  |i,  and  then  taking  its  derivative  with  respect  to  y,  one 
obtains  the  right-hand  side  of  the  momenmm  equation; 


(2.22) 


and 


puoo^  8 


24H  an 


pooiioo  ^yo j  an^ 


(2.23) 


The  nondimensional  form  of  the  momentum  equation  is  nearly  complete.  Setting  the  left- 
and  right-hand  sides  of  the  momentum  equation  equal  to  each  other,  one  obtains 


PUoo 


j_  a2f  ^  an  ^ 
a^an  ^  ^  an^. 


af  +  — — +  2&H^  — — 

an'2^HL  dn2  ^a^an2  ^  ^anan^ 


puoo^  a 

r  p 

r  >2n  a2f  ■ 

”  2^  an 

.PooMoo  [y 

oj  an^. 

^g8o(P<»'p) 


(2.24) 


Equation  (2.24)  can  be  simplified  by  multiplying  both  sides  by  (2^)/(puoo  ),  and  by 
expressing  the  first  and  second  derivatives  of  f  (with  respect  to  n)  as  f  and  f, 
respectively.  The  transformed  momentum  equation  becomes 
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B  fr+  2^aggoH(p^-p)  _  ,,  3F  3f  ^ 


f  ( 

an 

,  Poo  Jioo  V 

(2.25) 


The  buoyancy  term  on  the  left-hand  side  can  be  expressed  in  terms  of  Grashof  number, 
Gr,  and  Reynolds  numbers.  Re.  These  dimensionless  parameters  are  defined  as 


aggoH3(p«yp^-l) 


(2.26) 


(2.27) 


where  Voo  =  Poo/poo.  Substituting  the  expressions  for  Gr  and  Re  into  the  buoyancy  term, 
one  obtains 


2^aggoH(poo-p)  ^  (p«.  -  pwVpw  ,  H2  ^  Voo2  _  Gr  (Poo  -  p)/p 
pu«,^  (poo-pw)/Pw  Voo2  Re2(p^.p)/p^ 


(2.28) 


This  buoyancy  term  can  be  expressed  in  terms  of  the  Richardson  number,  which  is  the 
ratio  of  natural  to  forced  convection  (buoyancy  to  inertial  forces),  and  is  defined  as 


(2.29) 


The  final  form  of  the  transformed  momentum  equation  is 


CT"  f”  V  ff”  +  24Ri  ~  =  2  ^  ff  —  -  —  (2.30) 

anl^poopooW  J  (Poo-pw)/pw  V  a^  a^  ) 
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Mixture  Fraction  Equation 

In  a  similar  fashion,  the  transformed  equation  for  mixture  fraction.  Z.  can  be 
derived.  On  the  left-hand  side  of  Eq.  (2.3),  u  and  v  are  already  known.  Taking  the 
derivatives  of  Z  with  respect  to  x  and  y,  one  obtains 


8x  ■  ^  ax  ax  ^  a^  h 


and 


az_azan  ^  aza4_azr  uoo  im  ^ 
^  an  ^  a^  ^  an  L2^poomJ  ^  v" 

^  a 


(2.32) 


Upon  substituting  Eqs.  (2.31)  and  (2.32,  and  the  expressions  for  u  and  v,  the  left-hand 
side  of  Eq.  (2.3)  becomes 


az  az 


^an 

Lan^'" 


azj_ 

a^H 


oo  PooU  oo  ^/2 

n 


py 


V 


2^H 


f  +  2;— +  2^H  — ^ 


_ 1 

» 

leu 

In 

r  Uoo  *11 

J  ^Tl 

.2^1^  poo|-ioaj 

(2.33) 


or 


az  az 

p“ar-"p''a7= 


pUoo  „  3Z 


P^OQ  3Z  pUoo  3Z 

2^H  ^ 


(2.34) 
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Similarly,  the  right-hand  side  terms  of  Eq.  (2.3)  can  be  derived.  Multiplying  dTJdy  by  y", 
p,  and  D,  and  then  taking  the  derivative  with  respect  to  r\,  the  following  equation  is 
obtained: 


y"an 


ynpD|? 


n 

yo 


dri 

57 


(2.35) 


where  the  binary  diffusion  coefficient,  D,  is  expressed  in  terms  of  the  Schmidt  number. 
Sc,  such  that 


D  = 


(2.36) 


The  unity  Lewis  number  assumption  (Sc  =  Pr)  is  invoked  in  the  above  equation  to  express 
D  in  terms  of  the  Prandtl  number,  Pr.  Substituting  this  result  into  Eq.  (2.35)  and 
performing  some  manipulations,  the  right-hand  side  of  the  mixture  fraction  equation 
becomes 


pp  1  y"  dz 

[y  2^H  5ti 

pooPoo^yj  an 

(2.37) 


The  left-  and  right-hand  sides  of  the  mixture  fraction  equation  are  equated  to  obtain  the 
final  form  of  the  transformed  mixture  fraction  equation. 


3  r  P  ^  Wy  yn 
dr)  ^poo  |Aoo  ^  \y®J 


5nJ  an  Ha^a^anJ 


(2.38) 
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The  boundary  conditions  are  specified,  satisfying  the  non-slip  wall  and  uniform- 
flow  ambient  conditions: 

^>0  Ti=0  f=fw,f  =  0,Z=l 

-♦oo  f  =  l,Z  =  0.  (2.39) 

The  solution  at  ^  =  0  is  found  by  setdng  ^  to  zero  for  the  right-hand  sides  of  Eqs.  (2.25) 
and  (2.38).  Similarity  solutions  are  obtained  at  ^  =  0. 


Fuel  Burning  Rate 

An  expression  for  the  fuel  burning  rate  per  unit  surface  area  can  be  derived  from 
the  expression  (see  Lyu,  1991)  for  the  wall  velocity,  Vw, 


kw  1  aYpwazdTil 

'''“'p»cp.YF*-i“^a„37L 


(2.40) 


where 


is  evaluated  by 


an 

37 


(2.41) 


Substituting  the  above  equation  into  Eq.  (2.40),  Vw  becomes 


vw  = 


kw _ 1  avpw  az  Uoo 

PwCpw^^'^  ^  an  .2^Hpooiioo. 


(2.42) 
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The  termf — - ^1  is  multiplied  by  ^  and  —  to  yield  an  expression  that  includes  the 

L2^HpooMooJ  H  ^ 

Reynolds  number.  Re,  based  on  height  H.  The  equation  for  the  local  fuel  burning  rate 
then  becomes 


m"  pwVw  = 


1  pwkw  1  /Re^l/2^YFw^Z| 

^^2Hp.Cp*'’Fw-l[5  J  "5z 


(2.43) 


The  local  fuel  burning  rate  is  integrated  over  the  entire  wick  surface  to  obtain  the  total 
fuel  burning  rate,  given  by  the  following  equation: 


m  =  H(l) 


(2.44) 


where  the  surface  area  considered  is  the  product  of  the  height  H  and  unity  (1),  regardless 
of  geometry.  This  is  done  so  that  the  total  fuel  burning  rates  for  both  cylindrical  and 
planar  wicks  can  be  compared. 


Heat  Transfer 

In  addition  to  the  fuel  burning  rate,  an  important  parameter  to  heat  transfer  is  the 
Nusselt  number,  Nu.  The  local  Nusselt  number  is  defined  as 


Nux 


h  X 

kw 


(2.45) 


where  h  is  the  convective  heat  transfer  coefficient.  The  convective  heat  transfer 
coefficient  is  obtained  by  assigning  a  characteristic  temperature  difference  which,  in  the 
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present  study,  is  defined  to  be  TsfTw.  where  T'sj  is  the  stoichiometric  tempemture  ^nd 
is  the  wall  temperature.  The  heat  transfer  rate  is  obtained  by  multiplying  the  enthalpy  of 
fuel  vaporization,  hfg,  with  the  fuel  mass  burning  rate,  m.  The  heat  transfer  coefficient 
can  then  be  defined  as 


h  = 


q/A  rii"  hfg 
Tsi-Tw 


(2.46) 


where  m"  is  the  fuel  burning  rate  per  unit  surface  area.  The  local  Nusselt  number  is 


Nux 


m"  hfg  ^  H 
l^w  (TsfTw) 


(2.47) 


The  average  Nusselt  number,  Nujj,  over  the  surface  of  the  wick  is  defined  by 


hH 

kw 


(2.48) 


where  h  is  the  average  convective  heat  transfer  coefficient: 
h  =  ;^  JhdAs 

S  A, 


(2.49) 


In  general,  the  wick  surface,  Ag,  is  expressed  as 


(2.50) 
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where  n  is  0  for  a  flat  plate  and  1  for  a  cylinder.  The  derivative  of  As  is  taken  with 
respect  to  ^  to  obtain  an  expression  for  dAg: 


dAs  =  H(2jcyo)"d4 


(2.51) 


Equations  (2.46),  (2.50),  and  (2.51)  are  substituted  into  Eq.  (2.49): 


H(27ty, 


1 _  r 

Ttyo)"  J 


f&H(27cyoy'd^ 

1  yjj 


(2.52) 


Simplifying  Eq.  (2.52),  one  obtains 


‘'  =  H(Tst-Tw)  J“" 


(2.53) 


HCTst-Tw) 


(2.54) 


Finally,  substituting  Eq.  (2.54)  into  (2.48),  one  obtains  the  expression  for  the  average 
Nusselt  number 


NuH~if  V 

”  kw  (TsfTw) 


(2.55) 
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CHAPTERS 

NUMERICAL  SOLUTION  PROCEDURE 

This  chapter  describes  the  procedure  used  to  obtain  numerical  solutions  to  the 
governing  equations.  The  procedure  is  centered  around  a  computer  program  that  solves 
the  transformed  governing  equations.  The  code  requires  that  the  state  relationship  for  Li- 
SFg  combustion  be  specified.  The  computer  program  used  in  this  study  was  modified 
from  Lyu  (1991)  to  handle  mixed  convective  burning. 

To  obtain  solutions  to  the  governing  equations,  the  state  relationship,  e.g.  p  = 
p(Z),  must  be  specified.  The  state  relationship  was  obtained  from  the  equilibrium 
calculation  employing  the  NASA  code  with  thermodynamic  properties  of  Li2S  estimated 
from  Groff  (1976).  As  one  example,  the  state  relationships  for  temperature  and  density, 
taken  from  Lyu  (1991),  are  shown  in  Figure  3-1.  (The  reader  should  note  that  the  figures 
are  at  the  end  of  the  chapter.)  All  the  state  relationships  used  in  the  flow  calculations  are 
tabulated  in  Appendix  A,  c.f.  Table  A-1.  Note  that  the  interface  temperature  corresponds 
to  Z  =  1.  These  state  relationships  are  based  on  the  converged  interface  condition  (P  = 
0.01  MPa,  Ypw  =  0.38),  and  are  used  in  the  present  study  to  obtain  the  fuel  burning  rates 
at  0.01  MPa. 

The  flow  calculation  program  used  by  Lyu  (1991),  which  is  for  natural  convective 
burning,  was  modified  to  account  for  forced  convective  burning,  following  the  approach 
of  Chen  and  Faeth  (1981).  The  modifications  to  Lyu's  code  underwent  several  steps. 
First,  the  code  was  converted  to  pure  forced  convective  burning.  Gravity  was  no  longer 
considered  (i.e.,  ag  =  0).  Certain  coefficients  in  the  transformed  governing  equations 
were  changed,  as  well  as  the  boundary  conditions.  The  program  was  tested  for  a  flat 
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plate,  using  constant  properties  and  neglecting  wall  blowing  effects.  The  Blasius  solution 
was  obtained,  as  shown  in  Figure  3-2. 

The  code  was  then  modified  to  handle  mixed  convective  burning.  A  listing  of  the 
FORTRAN  computer  code  is  included  in  Appendix  A.  The  flow  calculation  program 
utilizes  a  finite  difference  scheme  (Keller  box  method)  in  solving  the  governing  equations 
at  specified  boundary  conditions  (Cebeci  and  Bradshaw,  1977).  A  marching  procedure  is 
used  to  obtain  the  solution  along  the  streamwise  distance.  The  computation  employs  82 
and  600  nodal  points  in  the  i)  and  ^  directions,  respectively.  Variable  grid  sizes  are  used 
for  both  directions.  The  first  100  ^  grid  points  are  spaced  1  x  10'^  apart,  and  then  spaced 
2  X  10"^  apart  for  the  remaining  500  grid  points.  The  T|  grid  points  are  separated  by  a 
value  Ati,  which  grows  with  each  successive  grid  generated;  each  successive  value  of  Aq 
is  multiplied  by  a  factor  of  1.05  as  the  grid  is  generated  from  the  surface  (i.e.,  qo  =  0,  qi 
=  1  X  10-2,  Aq2  =  1.05  •  1  x  10-2,  Aqs  =  1.05^  •  1  x  10-2,  ...  ,  Aqi  =  1.05*'^  •  1  x  10-2). 
This  allows  for  a  more  dense  grid  near  the  surface  and  a  less  dense  grid  farther  from  the 
surface. 

In  the  process  of  achieving  grid  independence,  the  nodal  mesh  required 
refinement  due  to  the  results  of  the  calculations  for  the  local  fuel  burning  rate;  the  results 
shows  some  peculiarities.  Specifically,  instead  of  decreasing  "uniformly"  along  the 
distance  of  the  surface,  the  local  fuel  burning  rate  showed  a  "stepping"  decrease  along  the 
streamwise  distance.  Figure  3-3  illustrates  the  stepping  phenomenon  in  more  detail. 
Further  investigation  of  the  stepping  change  found  that,  in  the  equation  for  the  local  fuel 
burning  rate,  the  gradient  of  the  mixture  fraction  was  fluctuating. 

The  stepping  change  was  reduced  when  the  grid  size  along  the  wick  surface  was 
refined.  Figure  3-4  shows  the  results  of  how  varying  the  grid  size  in  the  ^  direction 
affects  the  mixture  fraction  gradient  at  the  wall,  Z’lw,  which  consequently  affects  the 
value  for  the  local  fuel  burning  rate;  refer  to  Eq.  (2.43).  Four  grid  sizes  were 


22 


investigated,  corresponding  to  254,  504,  600,  and  1000  nodal  points  in  the  ^j-direction. 
The  first  case,  which  used  254  nodal  points  in  ^-direction,  is  the  same  grid  Lyu  (1991) 
used.  The  first  five  grid  points  were  variably  spaced,  while  the  remaining  249  grid  points 
were  spaced  4  x  10’3  apart,  beginning  at  ^  =  4  x  10‘3.  Since  the  stepping  changes  seemed 
more  pronounced  near  the  leading  edge  (^  <  1  x  lO’^)  and  since  the  boundary  layer 
thickness  changes  more  rapidly  near  the  leading  edge,  more  attention  was  given  to  this 
area  of  the  domain.  In  the  second  case  for  504  nodal  points,  sizes  of  the  first  four  grids 
(first  five  nodal  points)  were  reduced  to  one-fourth  the  size  of  those  in  the  first  case. 
Beginning  at  £  =  1  x  10'3,  the  remaining  500  nodal  points  were  spaced  uniformly  by  a 
distance  of  2  x  10'^.  The  third  case  involved  dividing  the  region  of  ^  <  1  x  10*^  into  99 
grids  of  size  1  x  lO'^  (100  nodal  points).  Finally,  in  the  founh  case,  1000  nodal  points 
were  employed:  500  nodal  points  spaced  2  x  lO'^  apart  (^  <  1  x  10*3)  and  500  nodal 
points  spaced  2  x  10*3  apart  (^  >  1  x  10*3). 

The  error  in  calculating  the  local  fuel  burning  rates  can  be  estimated  from  the 
results  for  Z'l^.  For  the  case  of  forced  convective  burning  (Re  =  3800)  along  a  flat  plate, 
the  fluctuation  ranges  from  a  minimum  of  -1.7  to  a  maximum  of  -1.85,  resulting  in  a 
maximum  difference  of  about  8%.  The  results  are  shown  in  Figure  3-5.  The  maximum 
fluctuation  occurred  immediately  after  ^  =  0.001,  where  the  grid  spacing  jumped  from  1  x 
10'^  to  2  X  10'^.  The  fluctuation  can  be  damped  out  by  refining  the  grid.  In  using  2000 
grid  points  with  variable  grid  spacing  (A^n  =  1  005  •  A^n-l).  Figure  3-6  shows  little 
variation  in  Z'  at  the  wall.  Indeed,  the  value  for  Z’  at  the  wall  andrh"  are  sensitive  to  the 
grid  scheme  used,  but  the  total  fuel  burning  rate,  th,  is  not  effected  after  grid 
independence  is  achieved.  The  calculated  total  fuel  burning  rates  for  the  cases  of  600  and 
2000  grid  points  differed  by  less  that  1%. 

The  total  fuel  burning  rates  were  obtained  by  performing  a  numerical  integration 
of  the  local  values  along  the  surface,  from  to  1,  using  Eq.  (2.44).  Simpson's  rule  was 
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used  for  the  integration  (Weltner  et  al.,  1986).  The  total  fuel  burning  rate  was  calculated 
for  different  values  of  ^min-  The  lower  integration  limit  was  chosen  to  be  10“^,  using  the 
results  shown  in  Figure  3-7  for  forced  convective  burning  along  a  flat  plate.  Between 
^min  =  10  ^  and  10^,  there  is  a  0.42%  change.  A  0.26%  change  occurs  between  10"^  and 
10-5 


aE 
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Mixture  Fraction 


Figure  3-1.  State  relationship  of  temperature  and  density  of  Li-SF6  wick  combustion 

Ypw  =  0.38,  taken  from  Lyu  (1991). 


Local  Fuel  Burning  Rate,  m"  (kg/sW  ) 
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Figure  3-2.  Blasius  solution. 


Dimensionless  Streamwise  Direction 


Figure  3-3.  Stepping  phenomena  illustrated:  local  fuel  burning  rate  along  planar  wick 
surface,  forced  convective  burning,  Re=3800  (Uoo=l  m/s,  H=1(X)  mm,  ag  =  0). 
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Figure  3-4.  Fluctuations  of  mixture  fraction  gradient  at  the  wall  for  various  grid 
schemes:  254,  504, 600, 1000  nod^  points  (^-direction). 
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Figure  3-5.  Semi-log  plot  of  fluctuation  of  mixture  fraction  gradient  at  the  wall,  forced 
convective  burning  along  a  flat  plate,  Re=3800;  600  nodal  points  (^-direction). 
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Figure  3-6.  Semi-log  plot  of  fluctuation  of  mixture  fraction  gradient  at  the  wall,  forced 
convective  burning  along  a  flat  plate,  Re=3800;  2000  nodal  points  (^-direction). 
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Figure  3-7.  Determination  of  minimum  integration  limit  from  total  fuel  burning  rate. 
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CHAPTER  4 

RESULTS  AND  DISCUSSION 

Though  presenting  the  results  for  mixed  convective  burning  is  the  main  objective 
of  this  study,  examining  forced  convective  burning  offers  valuable  insight  into  mixed 
convective  burning.  Through  a  parametric  study,  this  chapter  examines  forced  and  mixed 
convective  burning  of  Li-SF6  wick  diffusion  flames.  Both  wick  geometry  and  ambient 
conditions  are  examined  for  their  effects  on  burning  rates.  Cylindrical  wick  geometry  is 
emphasized,  since  it  yields  enhanced  burning  rates.  The  parameters  varied  in  the  study  of 
mixed  convective  burning  are  the  free-stream  velocity  of  the  SF6  oxidant  (Uoo),  the  wick 
height  (H),  the  gravity  (ag),  and  the  aspect  ratio  (heigh t/radius). 

Forced  Convective  Burning 

Mixed  convective  burning  combines  both  natural  and  forced  convection.  To 
develop  a  better  understanding  of  mixed  convection,  one  may  view  its  components 
separately,  investigating  the  burning  characteristics  associated  with  each  component.  In 
this  section,  the  results  for  forced  convective  burning  (ag  =  0)  are  presented  for  both 
planar  and  cylindrical  geometries.  In  the  study  of  forced  convective  burning,  the 
parameter  subject  to  variation  is  the  free-stream  velocity,  Uoo;  however,  the  results  are  in 
terms  of  the  Reynolds  number,  as  defined  by  Eq.  (2.27),  based  on  the  wick  height,  H. 
Other  parameters,  such  as  wick  height  and  gravity,  remain  constant  in  the  study.  The 
wick  height  is  100  mm,  and  ag  is  0.  For  cylindrical  wicks,  the  radius  is  fixed  at  12.5  mm. 
The  system  pressure  considered  in  the  computation  is  0.01  MPa.  The  state  relationships 
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and  interface  condition  were  obtained  from  Lyu  (1991).  The  interface  had  a  lithium  mass 
fraction  of  0.38  and  a  temperature  of  1340  K. 

Local  fuel  burning  rates  for  planar  and  cylindrical  wick  diffusion  flames  are 
calculated  by  Eq.  (2.43).  The  results  are  presented  in  Figures  4-1  and  4-2  for  Re  ranging 
from  3800  (Uoo  =  1  m/s)  to  190,000  (50  m/s).  A  dependence  is  observed  for  the 
planar  case.  The  results  for  the  cylindrical  case  deviate  from  the  dependence, 
however,  due  to  the  curvature  effects.  The  fuel  burning  rates  for  natural  convective 
burning  (flat  plate  and  cylinder)  under  the  same  ambient  and  geometry  conditions  (i.e.,  P 
=  0.01  MPa,  H  =  100  mm,  and  yo  =  12.5  mm)  are  also  graphed  for  comparison.  In 
contrast  to  forced  convective  burning,  natural  convective  burning  of  a  Li-SFe  wick 
diffusion  flame  shows  a  dependence  for  flat  plates  and  a  dependence  for 
cylinders  (Lyu,  1991), 

The  dependence  is  typical  for  the  heat  and  mass  transfer  across  laminar 
forced  convection  boundary  layers  (Williams,  1985).  For  wick  diffusion  flames,  the 
dependence  of  v^  is  a  sufficient  and  necessary  condition  for  the  right-hand  side  of  Eqs. 
(2.30)  and  (2,38)  to  vanish,  i.e.  fw  =  constant.  The  mathematical  derivations  are  included 
in  Appendix  B.  The  solution  yields  similar  profiles  for  f,  f,  f,  Z,  and  Z'  at  different 
streamwise  locations  when  the  similarity  variables  are  employed.  Figure  4-3  shows  a 
representative  sample  of  similarity  profiles  of  the  dimensionless  streamwise  velocity,  f . 

When  curvature  effects  are  considered,  an  additional  factor,  (y/yo)  ,  appears  in  the 
diffusional  transport  term.  The  local  fuel  burning  rates  of  cylinders  deviate  from  the 
dependence  of  flat  plates  for  %  >  0.001.  Similarity  solutions  are  not  obtained  in  the  case 
of  the  cylindrical  wick  due  to  the  curvature  term.  Overall,  the  fuel  burning  rates  for  the 
cylinder  are  slightly  higher  than  those  for  the  flat  plate.  According  to  Lyu  and  Chen 
(1991),  under  natural  convective  burning  conditions  the  curvature  effect  on  m"  becomes 
important  when  yo  «  H,  or  when  needle-type  cylinders  are  considered.  In  addition,  as  the 
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cylinder's  radius  is  increased  (keeping  the  height  constant)  m"  is  observed  to  approach  the 
flat  plate  solutions.  For  example,  when  the  cylinder  radius  yo  increased  to  50.8  mm  (H  = 
100  mm),  the  local  fuel  burning  rate  of  an  ethanol-air  wick  diffusion  flame  approaches 
that  of  a  planar  geometry  (Lyu  and  Chen,  1991).  Similar  behavior  is  observed  in  the 
present  study  of  Li-SFe  wick  diffusion  flames.  Using  the  same  wick  dimensions  as  Lyu 
and  Chen  (1991),  the  total  fuel  burning  rates  for  Li-SF6  cylindrical  wick  diffusion  flames 
under  forced  convection  conditions  approach  the  results  of  a  flat-plate  wick  geometry. 

The  total  fuel  burning  rates,  m,  were  calculated  using  the  corresponding  local  fuel 
burning  rates  calculated  for  various  Re.  This  was  done  by  performing  a  numerical 
integration  over  the  height  of  the  wick  surface.  As  discussed  earlier,  Simpson's  rule  was 
used  (Weltner  et  al.,  1986),  with  the  lower  integration  limit,  ^min.  set  at  0.0001.  The 
results  for  the  planar  and  cylindrical  cases  are  shown  in  Figure  4-4.  The  total  fuel 
burning  rates  for  the  cylindrical  flames  are  higher  than  those  for  the  flat  plate  flames;  m  is 
proportional  to  Re®'^  and  Re®'^^  for  the  flat  plate  and  cylinder,  respectively.  These 
relationships  are  true  for  the  parameters  specified.  The  total  burning  rate's  dependence  on 
Re  changes  if  a  different  wick  geometry  is  considered. 

Mixed  Convective  Burning 

The  results  from  the  previous  section  are  for  pure  forced  convective  burning,  i.e., 
ag  =  0;  buoyancy  effects  are  neglected,  and  Gr  is  set  to  zero.  In  reality,  however,  low 
gravity  can  be  expected  even  in  experimental  conditions,  such  as  in  a  space  station,  where 
typically  0  <  ag  <  0.01.  Furthermore,  during  the  technology  development  for  liquid-metal 
wick  burners,  buoyancy  effects  will  be  present  in  earth-bound  experiments.  At  these 
conditions,  mixed  convective  burning  needs  to  be  examined,  as  to  extrapolate  the  data 
base  for  space  applications. 
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The  following  section  examines  the  effects  of  changing  particular  parameters  on 
mixed  convective  burning.  These  parameters  are  the  free-stream  velocity  (Uoo),  the  wicK 
he’ght  (H),  the  arbitrary  gravitational  constant  (ag),  and  the  aspect  ratio  (H/yo).  The 
results  for  changing  the  free-stream  velocity  are  expressed  in  terms  of  Richardson 
number,  Ri,  instead  of  Uoo,  for  sake  of  clarity.  The  effects  on  the  burning  rates  and  heat 
transfer  for  both  cylindrical  and  planar  wick  geometries  will  be  highlighted. 


Free-stream  Velocity 

Mixed  convection  combines  both  natural  and  forced  convection.  The  important 
parameter  used  in  mixed  convection  is  the  Richardson  number  (Ri),  which  is  the  ratio  of 
natural  to  forced  convection.  As  stated  earlier  in  Eq.  (2.29),  Ri  is  defined  as 


Ri  can  be  stated  more  simply  by  considering  the  definitions  of  Gr  and  Re  and  then 
dividing  and  eliminating  terms.  The  following  expression  for  Ri  is  obtained: 

Ri  =  ^sSi£rPff.i]  (4  1) 

Uoo2  j 

From  Eq.  (4.1)  one  can  notice  that  Ri  decreases  as  Uoo  is  increased  by  a  power  of  2.  In 
this  computation,  the  other  parameters  are  fixed  (ag  =  0  and  H  =  100  mm). 

In  mixed  convective  burning,  forced  and  natural  convection  take  turns  dominating 
along  the  height  of  the  plate.  Near  the  leading  edge,  forced  convection  dominates. 
Figures  4-5  and  4-6  illustrate  this  phenomenon  using  the  dimensionless  streamwise 
velocity  profiles  (f  vs.  T|).  At  ^  =  0  of  a  flat  plate,  the  maximum  velocity  is  the  free- 
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stream  velocity  for  Ri  =  34  (uoo  =  1  m/s).  But  as  the  flow  moves  downstream,  buoyancy 
effects  come  into  play.  The  maximum  velocity  increases  as  ^  increases,  since  the 
buoyancy  forces  are  in  *he  same  direction  as  the  flow.  At  the  trailing  edge  the  peak 
velocity  is  almost  six  times  the  free-stream  velocity.  The  trends  are  similar  to  those 
obtained  for  aiding  mixed  convection  flow  over  a  vertical  surface  (Jaluria,  1986).  These 
trends  are  also  observed  in  the  case  of  the  cylinder.  The  peak  velocity  at  ^  =  1  is  over 
five  times  greater  than  the  free-stream  velocity. 

Burning  Rates 

Results  for  the  local  fuel  burning  rates  (m")  show  that,  as  Ri  approaches  zero 
(Uoo— >‘>®),  the  results  approach  those  for  forced  convective  burning.  Figures  4-7  and  4-8 
show  the  results  for  natural,  forced,  and  mixed  convective  burning  along  planar  and 
cylindrical  wicks,  respectively.  The  flat  plate  results  suggest  that  there  exists  a  critical  Ri 
beyond  which  m"  approaches  the  dependence  in  the  near-leading  edge  region,  or  ^  < 
0.01,  and  approaches  the  dependence  for  ^  t  0.01.  The  Grashof  number,  Gr, 
defined  by  Eq.  (2.26),  remains  the  same  (Gr  =  4.93  x  10^)  for  all  the  cases  involving 
mixed  convection.  The  total  fuel  burning  rates  for  the  burning  curves  of  Figures  4-7  and 
4-8  were  calculated,  and  the  results  are  shown  in  Figure  4-9.  Overall,  the  values  for  the 
cylinder  are  greater  than  those  for  the  flat  plate.  The  decrease  in  the  total  fuel  burning 
rate  (m)  as  Ri  increases  is  due  to  the  diminishing  contribution  of  forced  convection.  As 
mentioned  earlier,  Ri  increases  as  u«  decreases.  For  forced  convective  burning,  the  total 
fuel  burning  rate  decreases  as  Uoo  decreases.  Furthermore,  the  total  fuel  burning  rates  for 
mixed  convective  burning  are  higher  than  those  for  forced  convective  burning.  As  an 
example,  for  Uoo  =  10  m/s  (cylinder),  thmc  is  7.4  x  10"^  kg/s,  compared  to  6.9  x  10^  kg/s 
for  ihfc. 
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The  question  of  whether  the  results  for  forced  and  natural  convective  burning  can 
be  used  to  approximate  mixed  convective  burning  is  answered  by  the  results  tabulated  in 
Table  4-1.  As  the  previous  figures  showed,  ^or  the  region  near  the  leading  edge,  the 
mixed  convective  burning  curve  approaches  that  of  forced  convective  burning.  For  the 
region  closer  to  the  trailing  edge,  the  curve  approaches  that  of  natural  convective  burning. 
One  may  think  that  an  additive  relationship  among  the  different  modes  of  burning  might 
exist;  that  is,mnic  =  info  +  nine-  The  burning  rate  calculations  summarized  in  Table  4-1, 
however,  indicate  that  such  an  additive  relationship  does  not  exist 

Despite  the  fact  that  an  additive  relationship  does  not  exist  among  the  burning 
rates,  for  certain  "critical"  values  of  Ri,  mixed  convective  burning  can  be  approximated 
as  either  forced  or  natural.  These  critical  values  of  Ri  are  found  by  plotting  the  ratio  of 
fuel  burning  rates  (natural-to-mixed  and  forced-to-mixed)  versus  Ri  for  cylindrical  and 
planar  wicks.  This  is  done  in  Figures  4-10  and  4-11.  For  a  planar  wick  of  height  100 
mm,  Ricr  is  0.013  for  forced  convective  burning  and  155  for  natural  convective  burning. 
Likewise,  for  a  cylindrical  wick  of  height  1(X)  mm  and  radius  12.5  mm,  Ricr  is  0.028  for 
forced  convective  burning  and  160  for  natural  convective  burning.  The  critical  Ri  for 
natural  convective  burning  found  in  this  study  agrees  with  other  work  done.  Using  a  1% 
criterion  to  compare  Nu  for  natural  and  mixed  convection,  Yao  (1985)  determined  the 
value  for  Ricr.nc  to  be  500— the  same  order  of  magnitude  as  the  one  determined  in  this 
study— for  a  semi-infinite  vertical  flat  plate  undergoing  mixed  convection  heat  transfer. 
These  results  are  intended  to  assist  those  performing  earth-bound  experiments  related  to 
space  applications.  For  the  given  geometry  (H  =  100  mm  and  yo  =  12.5  mm),  the  total 
fuel  burning  rate  for  forced  convective  burning  can  be  approximated  using  the  results 
from  mixed  convective  burning  at  the  critical  Richardson  number,  Ricr.fc  with  an 
uncertainty  of  1%.  These  critical  Ri  change  when  considering  different  wick  sizes. 
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Heat  Transfer 

In  terms  of  heat  transfer  and  the  local  Nusselt  number,  Nux,  the  results  can  be 
represented  using  Eq.  (2.45).  The  value  of  Nux  is  based  on  the  difference  of  the 
temperatures  at  the  wall  (1340  K)  and  at  stoichiometric  conditions  (3120  K);  refer  to 
Figure  3-1.  The  results  are  presented  graphically  in  Figures  4-12  and  4-13.  Equation 
(2.45)  shows  that  Nux  is  proportional  to  ^  to  some  power,  which  depends  on  the 
relationship  between  the  local  fuel  burning  rate  and  For  forced  convective  burning,  in 
which  m"  is  proportional  to  Nux  is  proportional  to  Similarly  for  natural 
convective  burning,  in  which  m"  is  proportional  to  Nux  is  proportional  to  For 
the  case  of  the  flat  plate  (Figure  4-12),  such  relationships  can  be  seen.  For  the  curve  (Ri 
=  0.013,  forced  convection  dominant),  the  slope  is  approximately  1/2.  For  Ri  =  1 15, 
natural  convection  dominates  near  the  trailing  edge,  where  the  burning  curve  has  a  slope 
of  approximately  3/4.  This  observation  on  the  calculated  Nusselt  number  agrees  with  that 
of  Yao  (1985),  which  stated  that  in  mixed  convection  the  buoyancy  force  accelerates  the 
flow  near  the  surface  and  consequently  thins  the  boundary-layer,  thus,  a  higher  heat 
transfer  rate  is  obtained  at  downstream  locations.  The  average  Nusselt  numbers,  defined 
by  Eq.  (2.55)  are  calculated  and  the  results  are  shown  in  Figure  4-14.  Similar  to  the  total 
burning  rates,  the  average  Nusselt  numbers  increase  as  Ri  decreases. 

Already  being  investigated  as  a  heat  source  in  undersea  propulsion  systems 
(Hughes  et  al.,  1983,  and  Lyu  et  al.,  1990),  Li-SF6  wick  combustion  has  the  potential  of 
being  used  in  space  applications.  For  instance,  the  auxiliary  power  generation  system  for 
the  Space  Station  could  utilize  a  combustor  consisting  of  a  number  of  cylindrical  wicks 
undergoing  forced  convective  burning.  The  heat  output  could  be  controlled  by  either 
varying  the  free-stream  oxidant  velocity  or  changing  the  number  of  wicks  used.  For  the 
first  case  in  which  Uoo  is  varied,  a  relationship  between  the  heat  output  (per  wick)  and  Re 
needs  to  be  determined.  Such  a  relationship  has  been  determined  already  for  a  cylindric  1 
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wick  considered  in  this  study  (H  =  100  mm  and  yo  =  12.5  mm).  For  the  second  case,  in 
which  the  wick  dimensions  and  free-stream  velocity  are  fixed,  the  heat  output  of  a 
cylinde’^  '•'*ds  to  be  calculated.  Consider,  for  example,  one  cylindrical  wick  (0. 1  m  in 
height,  0.0125  m  in  radius)  under  forced  convective  burning  conditions  (P  =  0.01  MPa, 
Uoo  =  10  m/s).  The  predicted  total  fuel  burning  rate  for  this  case  is  5.4  x  lO'^  kg/s.  If  one 
kilogram  of  Li  is  assumed  to  yield  57,400  kJ  of  energy  (Lyu,  1991),  then  the  heat  output 
of  one  cylinder  is  0.31  kW.  A  power  system  requiring  a  5-kW  heat  source  would  need  17 
of  these  wicks  considered. 


Wick  Height 

It  would  be  of  interest  to  see  how  the  fuel  burning  rates  vary  with  changes  in  wick 
height.  The  previous  calculations  have  been  based  on  a  wick  height  of  100  mm.  When 
considering  mixed  convective  burning  along  wicks  of  different  heights,  one  must  keep  in 
mind  the  appropriate  length  scales.  Because  of  this,  the  results  for  the  local  fuel  burning 
rates  are  plotted  versus  the  local  Richardson  number,  Rix,  which  is  defined  as 


Ri,=2gi2Srp=., 

Uoo  j 


(4.2) 


instead  of  the  dimensionless  streamwise  distance,  Figures  4-15  and  4-16  show  this 
trend  for  local  fuel  burning  rates.  For  both  planar  and  cylindrical  wick  geometries,  the 
burning  curves  overlap  one  another,  unlike  the  previous  cases.  The  curve  for  the  shonest 
height  (H  =  0.01  m)  corresponds  to  the  lowest  range  of  local  Richardson  numbers-- 
roughly,  less  than  Rix  =  1.  The  curve  for  H  =  0.1  m  spans  across  the  highest  range  of  Rix. 
When  the  total  fuel  burning  rates  are  calculated  for  the  given  range  of  heights,  a 
relationship  is  observed  between  and  Ri.  For  the  flat  plate,  ih  is  proportional  to  Ri®‘^^ . 
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For  the  cylinder,  ihis  proportional  to  Ri  \  Figure  4-17  shows  the  total  fuel  burning 
rates.  These  relationships  are  limited  to  the  conditions  considered  (Uoo  =  1  m/s,  ag  =  1 ,  yo 
=  12.5  mm)  and  are  expected  to  change  for  different  conditions.  The  graph  indicates  that 
increasing  the  wick  height  (greater  Ri)  increases  the  total  fuel  burning  rate.  This  is  to  be 
expected  since  a  taller  wick  means  more  surface  area,  thus,  a  higher  burning  rate. 
Furthermore,  cylindrical  wicks  are  observed  to  yield  higher  total  burning  rates  over 
planar  wicks. 


Gravity 

The  effect  of  gravity  on  mixed  convective  burning  is  studied  by  varying  the  value 
of  ag.  The  gravity  conditions  are  changed  by  setting  ag  at  a  value  between  0  and  1;  i.e.,  0 
<  g  <  go.  where  go  is  the  gravity  at  sea  level.  The  burning  rate  calculations  are  performed 
for  ag  =  1,  0.5,  0.1,  0.01,  and  0.001.  Typical  gravity  in  an  environment  like  the  Space 
Station  would  have  an  ag  of  less  than  0.01.  In  the  calculations,  the  other  parameters,  H 
(100  mm)  and  Uoo  (1  m/s),  remain  fixed.  The  local  fuel  burning  rates  for  flat  plate  and 
cylinders  are  presented  in  Figure  4-18  and  4-19.  The  curves  for  natural  and  forced 
convective  burning  are  also  graphed  for  comparison.  Since  gravity  is  a  parameter  that 
affects  natural  convection  only,  the  local  fuel  burning  rates  are  approximately  the  same 
near  the  leading  edge  (^  <  0.001),  where  forced  convection  dominates.  Moving 
downstream,  the  curves  diverge;  the  curves  corresponding  to  smaller  values  of  ag 
approach  the  forced  convective  burning  limit  because  the  relative  influence  of  natural 
convection  diminishes.  The  buoyancy  forces  are  not  as  strong  because  of  the  decrease  in 
gravity. 

The  total  fuel  burning  rates  for  planar  and  cylindrical  Li-SFs  wick  combustion  are 
calculated  and  presented  in  Figure  4-20.  As  expected,  the  total  fuel  burning  rates  are 
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higher  for  cylindrical  wick  geometries  than  for  planar  wicks.  Overall,  for  mixed 
convective  burning  with  variable  gravity  conditions,  the  values  of  ih  approach  those  for 
forced  convective  burning.  At  ag  =  0.001,  the  total  fuel  burning  rate  is  1.46  x  10*^  kg/s 
(flat  plate)  and  3.41  x  lO'^  kg/s  (cylinder);  the  corresponding  total  fuel  burning  rates  for 
forced  convective  burning  are  98.4%  (flat  plate)  and  99.6%  (cylinder)  of  their 
corresponding  mixed  convective  burning  rates  at  ag  =  0.001  (Ri  =  0.034).  Thus,  the 
burning  rates  for  mixed  convective  burning  approach  the  forced  convective  burning  limit 
as  ag  is  decreased. 


Aspect  Ratio 

As  mentioned  earlier  in  the  section  dealing  with  forced  convective  burning,  the 
results  for  the  total  fuel  burning  rate  (cylinder)  approach  the  results  of  a  flat  plate  as  the 
radius  is  increased,  keeping  height  fixed.  To  determine  at  what  point  increasing  the 
radius  would  no  longer  affect  the  burning  rate,  the  ratio  of  the  total  burning  rates 
(cylinders  of  different  aspect  ratios  compared  to  the  flat  plate  solution)  are  plotted  versus 
aspect  ratio  in  Figure  4-21.  The  aspect  ratio  is  defined  as  the  wick  height  divided  by  the 
radius,  H/yo-  As  the  aspect  ratio  is  decreased  (increasing  radius),  the  total  burning  rate 
approaches  the  value  of  the  flat  plate  solution.  In  the  case  of  forced  convective  burning, 
an  aspect  latio  of  less  than  0.125  means  the  total  burning  rates  for  cylindrical  and  planar 
wicks  differ  by  1%.  For  mixed  convective  burning,  the  total  fuel  burning  rates  for 
cylindrical  and  planar  wicks  differ  by  1%  for  an  aspect  ratio  of  0.25.  Figure  4-22  shows 
the  local  fuel  burning  rates  for  cylindrical  wicks  of  radius  12.5,  50.8,  and  100  mm,  as 
well  as  the  corresponding  flat  plate  solution  (yo  — >  “»).  The  curve  for  yo  =  12.5  mm 
deviates  significantly  from  the  flat  plate  solution;  however,  the  curves  for  the  other  two 
cylindrical  cases  deviate  only  slightly  from  the  flat  plate  solution.  These  results  are 
similar  to  those  found  by  Lyu  and  Chen  (1991)  for  ethanol-air  wick  diffusion  flames. 


Local  tuel  Burning  Rate,  ih"  (kg/s-m  )  ^  Local  Fuel  Burning  Rate,  m"  (kg/s-ni^) 
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1.  Predicted  local  fuel  burning  rates  of  planar  Li-SF6  wick  diffusion  flame; 

H  =  100  mm. 
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Figure  4-2.  Predicted  fuel  burning  rates  of  cylinder  Li-SF6  wick  diffusion  flame; 

H  =  100  mm,  =  12.5  mm. 


Total  Fuel  Burning  Kate,  m  (kg/s) 


Dimensionless  Cross-stream  Distance 


Figure  4-3.  Similarity  profiles  of  dimensionless  streamwise  velocity  for  wick-type 
Li-SF6  diffusion  flame;  H  =  1(X)  mm,  Uoo  =  1  m/s. 


Re 


Figure  4-4.  Predicted  total  fuel  burning  rates  of  flat-plate  (H  =  1(X)  mm)  and  cylinder 
(H  =  100  mm,  yo  =  12.5  mm)  Li-SF6  wick  diffusion  flame  as  functions  of  Re. 


Dimensionless  Velocity,  r  ■f'  Dimensionless  Velocity,  T 


Dimensionless  Cross-stream  Distance 


.  Dimensionless  streamwise  velocity  of  flat-plate  Li-SF6  wick  diffusion  flame; 
Ri  =  34  (Uoo  =  1  m/s),  H  =  100  mm,  ag  =  1 . 


Dimensionless  Cross-stream  Distance 


Figure  4-6.  Dimensionless  streamwise  velocity  of  cylinder  Li-SF6  wick  diffusion  flame, 
Ri  =  34  (uoo  =  I  m/s),  H  =  100  mm,  yo  =  12.5  mm,  ag  =  1 . 


l.ocal  I'uel  Burnin}>  Kate,  m  (kgis-m  )  ^  Local  Fuel  Burning  Rale,  m"  (kgjs-m^) 
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Dimensionless  Streamwise  Distance 

7.  FYedicted  local  fuel  burning  rates  for  flat-plate  Li-SFe  wick  diffusion 
flame;  H  =  100  mm,  ag  =  1.0,  Gr  =  4.93  x  10*. 
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Figure  4-8.  Predicted  local  fuel  burning  rates  for  cylindrical  Li-SF^  wick  diffusion 
flame;  H  =  100  mm,  ag  =  1.0,  Gr  =  4.93  x  10*. 
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Figure  4-9.  Predicted  total  fuel  burning  rates  for  flat-plate  (H  =  100  mm)  and  cylinder 
(H  =  100  mm,  yo  =  12.5  mm)  Li-SF6  wick  diffusion  flame  as  functions  of  Ri; 

ag  =  1.0,  Gr  =  4.93  x  lO®. 


Table  4-1.  Comparison  of  total  fuel  burning  rates  (kg/s)  for  cylinder 
(  H  =  100  mm,  yo  =  12.5  mm). 


Ri 

nine 

riifc 

riinc  '*■  mfc 

mQjc 

%  Difference 

34 

0.00035 

0.00014 

0.00049 

0.00037 

32 

2.1 

0.00035 

0.00029 

0.00064 

0.00042 

52 

0.36 

0.00035 

0.00045 

0.00080 

0.00052 
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0.037 

0.00035 

0.00079 

0.0011 

0.00081 
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0.0014 
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Log  Ri 

Figure  4-10.  Critical  Richardson  numbers  for  natural  and  forced  convective  burning 
limits  (99%  accuracy);  flat-plate  Li-SF6  wick  diffusion  flame;  H  =  100  mm,  ag  =  1.0. 

Ricr.fc  =  0.013  and  Ricrjic  =  155. 


I.og  Ri 

Figure  4-11.  Critical  Richardson  numbers  for  natural  and  forced  convective  burning 
limits  (99%  accuracy);  cylindrical  Li-SF^  wick  diffusion  tlame; 

H  =  100  mm,  =  12.5  mm,  ag  =  1.0.  Rkr.ic  =  0.028  and  Ricrjic  =  160. 


Local  Nusselt  Number,  Nu^  ^  Local  Nusselt  Number,  Nu 
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12.  Predicted  local  Nusselt  number  for  flat-plate  Li-SF6  wick  diffusion  flame; 

H  =  100  mm,  ag  =  1.0. 
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Figure  4- 1 3.  Predicted  local  Nusselt  number  for  cylinder  Li-SF6  wick  diffusion  flame; 

H  =  100  mm,  yo  =  12.5  mm,  ag  =  1.0. 
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Figure  4-14.  Predicted  average  Nusselt  numbers  for  flat-plate  (H  =  100  mm)  and 
cylinder  (H  =  100  mm,  yo  =  12.5  mm)  Li-SF6  wick  diffusion  flame  as  functions  of  Ri; 

ag  =  1.0,  Gr  =  4.93  x  10*. 
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Figure  4-15.  Effect  of  varying  wick  height  on  local  fuel  burning  rate  for  planar  wick 

(H  =  0. 1 , 0.05,  0.0 1 );  Uoo  =  1  m/s,  ag  =  1 . 


Local  Fuel  Burning  Kate,  m"  (kgls-m^) 
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Figure  4- 16.  Effect  of  varying  wick  height  on  local  fuel  burning  rate  for  cylindrical 
wick  (H  =  0.1, 0.05, 0.Ol);  yo  =  12.5  mm,  Uoo  =  1  m/s,  ag  =  1. 


Figure  4- 1 7.  Predicted  total  fuel  burning  rate  of  Li-SFs  wick  combustion  as  function  of 
Ri  (H  =  0. 1,  0.075, 0.05, 0.025, 0.01  m);  yo  =  12.5  mm,  Um  =  Im/s,  ag  =  1.0. 
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Figure  4-18.  Effect  of  varying  arbitrary  gravitational  constant,  a®,  on  local  fuel  burning 
rate  for  planar  wick  (ag  =  1, 0.5, 0.1, 0.01, 0.001);  H  =  100  mm,  Uoo  =  1  mys. 
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Figure  4- 19.  Effect  of  varying  arbitrary  gravitational  constant,  ao.  on  local  fuel  burning 
rate  for  cylindrical  wick  (ag  =  1, 0.5, 0. 1, 0.01,  OTOl ); 

H  =  100  mm,  yo  =  lz.5  mm.  ito  =  1  m/s. 


Burning  Kate,  m  (kg/s) 

m  (Mat  Plate) 
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Ri 


-20.  Predicted  total  fuel  burning  rates  of  planar  and  cylinder  Li-SF6  wick 
diffusion  flames  as  function  of  Ri  (ag  =  1, 0.5, 0.1, 0.01,  0.001); 

H  =  100  mm,  yo  =  12.5  mm,  u*  =  Im/s. 


Figure  4-2 1 .  Ratio  of  total  fuel  burning  rates  (cylinder  to  flat  plate)  versus  (H/yo); 
variable  radius.  Uoo  =  1  m/s,  H=  100  mm,  a®  =  1;  forced  and  mixed  convective  burning  of 

Li-SF6  wick  diffusion  flames. 
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CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 

The  objective  of  this  study  is  to  extend  the  methodology  (developed  for  Li-SFe 
wick  diffusion  flames  under  natural  convective  burning  conditions)  to  investigate  forced 
and  mixed  convective  burning  of  Li-SFg  wick  diffusion  flames.  A  conserved  scalar 
approach  is  employed  to  model  the  wick  diffusion  flames,  and  the  interface  condition  is 
determined  by  the  conservation  of  mass,  species,  and  energy  at  the  interface. 
Computation  is  first  done  for  forced  convective  burning  of  Li-SF6  wick  diffusion  flames 
and  then  extended  to  mixed  convective  burning.  A  limitation  of  the  present  analysis 
comes  from  the  employment  of  boundary-layer  approximations.  The  analysis,  in  essence, 
considers  a  semi-inflnitely  long  wick.  The  end  effects  of  a  finite  height  wick  are  not 
considered. 

In  the  process  of  obtaining  numerical  results,  the  computer  code  from  Lyu  (1991) 
was  modified  to  calculate  the  local  fuel  burning  rate.  The  code  utilized  the  state 
relationships  based  on  the  converged  interface  condition  (P  =  0.01  MPa,  Ypw  =  0.38), 
which  were  already  calculated  by  Lyu  (1991),  using  the  NASA  CEC  equilibrium  code 
(Gordon  and  McBride,  1976).  The  computer  program  used  a  finite-difference  scheme 
and  marching  procedure  to  solve  the  transformed  governing  equations.  The  numerical 
grid  employed  600  points  in  the  ^-direction  and  82  points  in  the  ti -direction.  The  total 
fuel  burning  rates  were  calculated  using  Simpson's  rule  of  numerical  integration  over  the 
length  of  the  wick  surface  (from  ^min  =  0.0001  to  ^  =  1). 

Flow  calculations  involving  forced  convective  burning  were  performed  for  both 
planar  and  cyl  -idrical  wick  geometries.  Similarity  solutions  were  obtained  for  flat-plate 
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Li-SFg  wick  diffusion  flames,  i.e.,ih"  «  Because  of  curvature  effects,  however, 

nonsimilar  solutions  were  obtained  for  cylinder  Li-SF6  wick  diffusion  flames.  In  general, 
cylindrical  wicks  enhanced  the  burning  such  that  the  total  fuel  burning  rates  were  higher. 
Flow  calculadons  were  then  done  for  mixed  convective  burning.  The  results  for  the  local 
fuel  burning  rates  showed  that  forced  convection  dominates  near  the  leading  edge  and 
natural  convection  over  the  rest  of  the  surface  The  results  indicated  that  natural 
convection  has  a  large  contribution  to  mixed  convective  burning. 

A  parametric  study  is  done  on  the  effects  of  varying  the  free-stream  velocity  (Uoo), 
the  wick  height  (H),  and  the  gravity  (ag)  on  the  fuel  burning  rate.  The  starting  values  for 
each  parameter  were  1  m/s  for  Uoo,  100  mm  for  H,  and  1.0  for  ag.  Each  calculation  was 
done  for  both  planar  and  cylindrical  wick  geometries.  The  effect  of  increasing  Uoo  was  to 
increase  the  local  fuel  burning  rates.  For  the  planar  case,  as  Uw  was  increased,  the 
burning  curves  were  observed  to  have  forced  convective  burning  characteristics,  i.e.,  ih" 
Upon  further  investigation,  critical  Richardson  numbers  were  found  for  a  given 
geometry  and  gravity  condition  (H  =  100  mm,  yo  =  12.5  mm,  and  ag  =  1.0);  these 
numbers  indicate  when  mixed  convective  burning  can  be  approximated  by  either  forced 
and  natural  convective  burning.  When  the  wick  height  H  was  decreased  from  100  mm  to 
10  mm,  the  total  burning  rates  were  observed  to  decrease.  A  relationship  between  the 
total  fuel  burning  rates  and  log  Ri  was  obtained  for  both  flat  plates  and  cylinders.  The 
local  burning  curves  for  different  heights  overlap  one  another  when  the  burning  rate  is 
plotted  versus  the  local  Richardson  number.  When  the  gravity  was  decreased  from  1 .0  to 
0.001,  burning  rates  also  decreased.  The  values  of  the  local  fuel  burning  rates  near  the 
leading  edge,  where  forced  convective  burning  dominates,  remained  relatively 
unaffected;  however,  downstream,  where  natural  convection  dominates,  the  local  burning 
rates  decreased  noticeably  due  to  the  lower  gravity,  which  resulted  in  smaller  buoyancy 
forces.  As  ag  was  decreased,  the  burning  rates  approached  the  forced  convective  burning 
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limit.  Regarding  the  aspect  ratio,  the  fuel  burning  rates  for  cylindrical  wicks  (both  forced 
and  mixed  convective  burning)  approached  the  results  for  planar  geometry.  In  the  case  of 
mixed  convective  burning,  the  local  fuel  burning  rate  for  a  cylindrical  wick  of  height  100 
mm  and  radius  50.8  mm  ^eviated  only  slightly  from  the  flat  plate  solution. 

This  study,  though  not  a  comprehensive  treatment  of  mixed  convective  burning  of 
Li-SF6  wick  diffusion  flames,  has  accomplished  its  objectives.  For  one,  a  mathematical 
model  was  developed  for  wick-type  Li-SFe  diffusion  flames  undergoing  mixed 
convective  burning.  Secondly,  numerical  solutions  were  obtained  to  assess  the  effects 
due  to  geometry  (cylindrical  versus  planar  wicks,  varying  wick  heights  and  aspect  ratios) 
and  ambient  conditions  (free-stream  velocity  of  oxidant,  gravity)  on  the  wick  combustion 
of  Li  and  SFe.  Thirdly,  results  were  presented  that  will  add  to  the  knowledge  base  of 
earth-bound  combustion  experiments  performed  for  variable  gravity  conditions.  These 
results  include  establishing  a  relationship  between  the  burning  rate  and  the  free-stream 
velocity,  investigating  the  effects  of  reduced  gravity  on  the  wick  combustion,  and 
performing  a  sample  calculation  to  determine  the  number  of  wicks  needed  to  meet  a  heat 
output  requirement. 

Further  analysis  and  experimentation  are  needed  to  evaluate  the  uncertainties  of 
the  analysis,  as  well  as  confirm  the  results  predicted.  Suggested  future  work  includes 
extending  the  analysis  to  turbulent  flow,  incorporating  radiative  heat  transfer  into  the 
analysis,  analyzing  the  end  effects  of  a  finite  wick,  and  considering  the  effects  of  the 
presence  of  non-condensable  gases,  such  as  argon. 
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APPENDIX  A 

COMPUTER  PROGRAM  FOR  FLOW  CALCULATION 


Main  Program 


C  FILENAME:  FLOW.MC.C. V.U  1  (U_  DENOTES  FREE-STRE AM  VELOCITY) 

C  COMMENTS:  TfflS  IS  THE  FLOW  CALCULATION  PROGRAM  FOR  MIXED  CONVECTION 
C  •  CYLINDER:  CURV AT  NOT  EQUAL  TO  0  DIMENSIONS:  H=0. 1  M  AND  R= 1 2.5  MM 

C  *  VARIABLE  PROPERTIES:  LI  AND  SF6  SYSTEM  PRESSURE  =  0. 1  MPA 

C  *  YFW  =  0.38  (CONVERGED  INTERFACE  CONDITION) 

C  *  OUTPUT:  LOCAL  FUEL  BURi'nNG  RATE  RMBR 


k *»*»»•» »»***»***» 


IMPLICIT  REAL*8(A-H,0-Z)JNTEGER*4(I-N) 

COMMON  /INPTI/  A1,A2,A1Z,ALFA0ALFA1,IT,IFLAG,GR3UA. 

1  HPLATEJIADIUS,RKVISI.RE.UINFJU 

COMMON  /EOS2/  ZMF(IOI)  JIHO(101).TEP(101)31HOM(200),TEMP(200), 

1  DVC(101)J>R(10l)T)VCM(200).PRM(200)JC 

COMMON  /ETA3/  ETAE.NP6J4P.DETA(200)^TA(200)  A(200);'IPC,NNPJ4XT. 
1  NX,X(600).CVW.VW(2).VWT 

COMMON  /REST4/  F(200,2),U(200,2),V(200,2)2(200^)4>(200^). 

1  AM{200,2),B(200,2).BZ(200,2),BT(200) 

COMMON  /SOLV6/  DELF(2(X)),DELU(200)T>ELV(200),DELZ(200). 

1  DELP(200).ABSERR 


ITMAX  =  900 


C — 

C 

C 


NX  =  I 


SUBROUTINE  INPUT  SUPPLIES  FLOW  GEOMETRY  AND  PROPERTY  CONSTANTS 

AIA2:  MOMENTUM  EQUATION  PARAMETERS 

A17:  MIXTURE  FRACTION  PARAMETERS 

/»J.FA0ALFAI:  BOUNDARY  CONDITIONS  PARAMETER 

IFLAG:  CONTROL  PARAMETER 

GR:  GRASHOF  NUMBER 

HPLATE:  LONGITUDINAL  LENGTH 

RADIUS:  RADIUS  OF  CYLINDER 

RKVISI:  KINEMATIC  VISCOSITY  AT  INFINITY 

RJA:  MODIFIED  JACOB  NUMBER 

RE:  REYNOLDS  NU'MBER 

UINF:  VELOCITY  AT  INFINITY 

RI:  RICHARDSON  NUMBER  (GR/RE''2) 


CALL  INPUT 


SUBROUTINE  GRID  SET  UP  GRID  SYSTEM 

NXT:  TOTAL  GRID  POINT  IN  X-DIRECTION 


n  o  n  o 
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C  X(NX):  GRID  POSITION  IN  X-DIRECTION 

C  ETAE:  BOUNDARY  LAYER  THICKNESS 

C  ETA(J):  GRID  POSITION  IN  Y-DIRECTION 

C  DETA(J):  INCREMENT  OF  GRID  POSITION  IN  Y-DIRECTION 

C  VGP;  RATIO  OF  ADJACENT  GRIDS  IN  Y-DIRECTION 

C  NP;  TOTAL  GRID  POINT  IN  Y-DIRECnON 

CALL  GRID 

C  SUBROUTINE  IVZPL  PRODUCES  INITIAL  GUESS  PROHLES  AND  BOUNDARY 

C  CONDITIONS  FOR  FLOW  CALCULATION 

C  F(J.2):  DIMENSIONLESS  STREAM  FUNCTION 

C  U(J,2):  DIMENSIONLESS  STREAMWISE  VELOCITY 

C  V(J,2):  DERIVATIVE  OF  U(J) 

C  Z(J^):  MIXTURE  FRACTION 

C  P(J.2):  DERIVATIVE  OF  Z(J) 

CALL  IVZPL 
IFLAG=  I 
10  IT  =  0 

20  IT  =  IT-H 

IF(  IT  .LT.  ITMAX )  THEN 

C  SUBROUTINE  INUPD  INITIALIZE  AND  UPDATE  PARAMETERS  AND  CALCULATE 
C  CURVATURE  TERM 

C  AM(J,2):  MOMENTUM  SOURCE  TERM 

C  B(j;2):  PROPERTY  TERM  IN  MOMENTUM  EQUATION 

C  BZ(J,2);  PROPERTY  TERM  EN  MIXTURE  FRACTION  EQUATION 

CALL  INUPD 
CALL  COEF 
CALL  SOLV5 

IF(  ABSERR  .GT.  LD-06)  THEN 
GO  TO  20 

ELSEIF(  DABS(  V(NP.2) )  .GT.  LD-03  )  THEN 
CALL  GROW  I'H 
C  WRITE(6.1000)  NNP 

1000  FORMAT(1HO,  'ETAE  GROW  ',13,’  -POINTS  ADDED’) 

GO  TO  10 

ENDDF 

ELSE 

C  WRITE(6,2000)  NX,X(NX),DELV(1),DELP(1) 

2000  FORMAT(1HO,4HNX  =,I3,5X,3HX  =,E10.4,5X,25HITERATIONS  EXCEEDED 
1  ITMAX,5X,7HDELVW  =E10.4,5X,7HDELZW  =£10.4) 

ENDIF 


CORRECT  WALL  BLOWING  VELOCITY  ACCORDING  TO  ENERGY  AND  SPECIES 
WALL  BOUNDARY  CONDITIONS 


IF(  NX  .NE.  1)  THEN 
VWT  =  'VW(2) 

C  VW(2)  =  -(  P(l,2)  +  ZPWT  )/2.D-t-00*(  GR/X(NX)  )**2.5D-01/RJA 
VW(2)  =  -P(l,2)  *  (RE/X(NX))*»5.0D-01/RJA 

C  VW(2)  =  -(  P(l,2)  -t-  ZPWT  )/2.D-(-00*DSQRT(  RE/X(NX)  )/RJA 
ZPWT  =  P(  1,2) 


noonn  oooo 
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IF(  DABS(  VW(2)  -  VWT )  .GT.  l.D-08 )  THEN 
F(l,2)  =  F(l,l)*(  X(NX-1)/X(NX)  )**5.0D-01 

1  +  CVW*(  VW(1)  +  VW(2)  )/2.D■^W*(  X(NX)  -  X(NX  -  1) ) 

2  /X(NX)**5.0D-01 
GO  TO  20 

ELSE 

VW(1)  =  VW(2) 

ENDIF 

ENDIF 

CALL  OUTPUT 

IF(  NX  .LE.  NXT )  THEN 

GO  TO  10 

ENDIF- 

WRITE(6,514) 

514  FORMAT(1HO,'MAIN:  PROGRAM  HAS  EXECUTED’) 

STOP 

END 

SUBROUTINE  INPUT 

IMPLICIT  REAL’*8(A-H,0-Z),INTEGER*4(I-N) 

COMMON  /INPTl/  A1.A2A1ZALFA0ALFAUT,IFLAG,GRJRJA, 

1  HPLATEJIADIUSJIKVISLRE,UINFJ^I 

COMMON  /EOS2/  ZMF(101)  jlHO(101).TEP(101)  JIHOM(200),TEMP(200). 
1  DVC(101),PR(101).DVCM(200),PRM(200)JC 


Al  =  l.D+00 
A2  =  O.D+00 
A1Z=  l.D+00 

ALFAO  =  1 ,  ALFA  I  =  0  SPECIFY  FUNCTION  VALUE 
ALFAO  =  0.  ALFAl  =  1  SPECIFY  FUNCTION  DERIVATIVE 

ALFAO  =  1 
ALFAl  =  0 


READ  IN  CONSERVED  SCALAR.  TEMPERATURE.  DENSITY.  VISCOSITY,  AND 
PRANDTL  NUMBER  FROM  CEC72  CALCULATION  (FROM  BLOCK  DATA  AT  END) 

**  SUPPRESS  PROPERTIES  FOR  BLASIUS  SOLUTION  ** 


1000 

C 

C 

C 

C 

C20 

C — 

C 

C 

C 

C 

C 

C . 


READ(5.1000)(ZMF(J),TEP(J)J^HO(J).DVC(J),PR(J)J  =  1,101  ) 

FORMAT(F7.4.5X  JT.  1 ,5X£  1 1 .4.5X.E1 1 .4.5X,F6.4) 

DO20J=l,101 

RHO(J)  =  5.9696D-01 

DVC(J)=  L5710D-05 

PR(J)  =  0.7368D+00 

CONTINUE 


FLOW  AND  GEOMETRY  CONSTANTS 

YFW:  FUEL  MASS  FRACTION  AT  WALL 

CPW:  SPECIFIC  HEAT  AT  WALL 

RKMW:  THERMAL  CONDUCTIVITY  OF  MIXTURE  AT  WALL 

DYFDZ:  DERIVATIVE  OF  FUEL  MASS  FRACTION  W.R.T.  MIXTURE  FRACTION 


HPLATE  =  l.D-01 
RKVISI=  1.55725D-05/RHO(1) 

GR  =  9.80665D+00'»DABS(  RHO(1)/RHO(10I)  -  I.D+00  )* HPLATE **3.0+00 
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1  /RKVISVRKVISI 

YFW  =  3.84741D-01 
CPW  =  2.7935D+03 
RKMW  =  7.315D-02 
DYFDZ  =  7.6273D-01 

C  RJA  =  DSQRT(  2.D+00  )*HPLATE*(  l.D+00  -  YFW  )*RH0(1)*CPW/RKMW 

RJA  =  DSQRT(  2.D400  )*HPLATE*(  l.D+OO  -  YFW  )*RHO(l)*CPW/RKMW/DYFDZ 
RADIUS  =  1.25D-02 
UINF=  l.D+00 

RE=UINF*HPLATE/RKVISI 

RI=GRyRE/RE 

RETURN 

END 

SUBROUTINE  GRID 

IMPLICIT  REAL*8(A-H.O-Z).INTEGER*4(I-N) 

COMMON  /ETA3/  ETAEJsfP6  J^P  J)ETA(200)£TA(200)  A(200)  J'lPCJWP.NXT, 

1  NX,X(600),CVW,VW(2),VWT 

0====^==:============================ 

NXT  =  600 
DXP  =  2.D-03 
DXPO=  l.OD-05 
L=  101 

DETA(1)=  l.D-02 
VGP=  1.05D+00 
ETAE  =  l.OD+o. 

C  GRID  GENERATION  IN  X-DIRECTION 

^ _ _ _ _ _ _ _ _ _ _ _ _ 

X(l)  =  0.D+00 
DO  95  I=2X 
X(I)  =  X(I-l)  +  DXPO 
95  CONTINUE 

DO  100I  =  L+1JVXT 
X(I)  =  X(I-1)  +  DXP 
100  CONTINUE 

C - - - 

C  GRID  GENERATION  IN  Y-DIRECTION 

C - 

IF(  (VGP  -  l.D+00)  .GT.  l.D-03  ) THEN 

NP  =  DLOG(  ETAE/DETA(1)’*'(VGP  -  l.D+00)  +  l.D+00  ) 

1  /DLOG(VGP)  +  l.OOOlD+00 

NP6  =  DLOG(  6.D+00/DETA(l)*(VGP  -  l.D+00)  +  l.D+00 ) 

1  /DLOG(VGP)  +  1 .0001 D+00 

ELSE 

NP  =  ETAE/DETA(1)  +  1.0001  D+00 
ENDIF 

IF(NP.LE.  310)  THEN 
ETA(l)  =  O.D+00 
D0  200J  =  2.NP 
DETA(J)  =  VGP*DETA(J-1) 

A(J)  =  DETA(J-l)/2.D+00 
ETA(J)  =  ETA(J-l)  +  DETA(J) 

200  CONTINUE 
ELSE 

WRITE(6.1000) 
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1000  FORMATdHO.NP  EXCEEDED  MAXIMUM  ARRAY  DIMENSION  ALLOWED 
1  PROGRAM  TERMINATED') 

STOP 

ENDIF 

RETURN 

END 

SUBROUTINE  IVZPL 

IMPLICIT  REAL*8(A-H.O-Z)JNTEGER*4{I-N) 

COMMON  /INPTl/  A1,A2A12ALFA0ALFAUT.IFLAG,GRJUA, 

1  HPLATEJIADIUS.RKVISI.RE.UINFJU 

1  NX,X(600),CVW.VW(2).VWT 

COMMON  /REST4/  F(200.2),U(200,2).V(200,2)^(200,2)JP(200,2), 

1  AM(200,2),B(200,2).BZ(200,2),BT(200) 

COMMON  /SOLV6/  DELF(200),DELU(200),DELV(200),DELZ(200). 

1  DELP(200),ABSERR 

0=============================== 

ETA6  =  6.D+00  +  5.D-01 

C  GENERATE  INITIAL  PROFILE  BY  SOLVING  THE  INCOMPRESSIBLE  FLOW 

DO  100J=  1J^6 
ETAR  =  ETA(J)/ETA6 
ETAR2  =  ETAR-^ETAR 

F(J,2)  =  ETA6/4.D+00*ETAR2*(  3.D+00  -  ETAR2/2.D+00 ) 

U(J,2)  =  ETAR*(  L5D+00  -  ETAR2/2.D+00  ) 

V(J,2)  =  L5D+00*(  l.D+00  -  ETAR2  )/ETA6 
Z(J^)  =  LD+00  -  ETAR 
P(J.2)  =  -LD+00/ETA6 
100  CONTINUE 

C-— - — - - - - — - - 

C  INITIAL  VALUE  OF  BOUNDARY  CONDITIONS 

U(NP6^)  =  l.D+00 
Z(NP6.2)  =  O.D+00 
VW(1)  =  0.I>KX) 

VW(2)  =  O.D+00 

C - 

C  GENERATE  BETTER  PRORLES 

C  EFLAG  =  0.  ETAE  =  ETA6  =  6.D+00 
C  IFLAG  =1,  ETAE  =  ETAE 

C - 

IFLAG  =  0 
IT  =  0 

10  IT  =  IT+1 

C  INITIALIZE  AND  UPDATE  PROFILES 

CALLINUPD 

IF(  IT  .GT.  100)  THEN 

WRnE(6.1000) 

1000  FORMAT(1HO,'INCOMPRESSIBLE  FLOW  SOLUTION  DID  NOT  CONVERGE') 
STOP 
ELSE 

CALL  COEF 
CALL  SOLV5 
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IF(  ABSERR  .GT.  l.D-06 )  GO  TO  10 
ENDIF 

C - 

C  SET  UP  DUMMY  NODES 

D0  200J  =  NP6+I^ 

F(J,2)  =  F(NP6^) 

U(J.2)  =  U(NP6^) 

V(J.2)  =  V(NP6^) 

Z(J^)  =  Z(NP6.2) 

P(J,2)  =  P(NP6^) 

200  CONTINUE 

RETURN 
END 

SUBROUTINE  COEF 

IMPLICn'  REAL*8(A-H.O-Z)  JNTEGER*4(I-N) 

COMMON  /INPTl/  A1,A2^1ZALFA0ALFA1.IT,IFLAG,GRJUA, 

1  HPLATEJIADIUSJIKVISITIE.UINFJU 

COMMON  /ETA3/  ETAEJ4P6J4PT)ETA(200)£TA(200)  A(200);'IPC.NNPJ«:T, 
1  NX,X(600).CVW.VW(2),VWT 

COMMON  /REST4/  F(200,2),U(200,2),V(200,2)^(200.2)  J’(200,2), 

I  AM(200,2)3(200,2),BZ(200,2),BT(200) 

COMMON  /COEF5/  S  1(200).S2(200).S3(200),S4(200).S5(200),S6(200), 

1  S7.S8.B  1(200), B2(200).B3(200).B4(200), 

2  B5(200)36(200).B7(200).B8(200)3l(5,200) 
0============================ 

IF(  IT  £Q.  1)  THEN 
CEL  =  0.D+00 
IF(  NX  .GT.  1)  THEN 

CEL  =  I.D+00*(  X(NX)  +  X(NX  - 1) )/(  X(NX)  -  X(NX  -  1) ) 

ENDBF 

P1P  =  AI  +  CEL 
P2P  5  A2  +  CEL 
P1PZ  =  AIZ  +  CEL 
ENDIF 

C - - - 

C  PRESENT  STATION 

C - 

DO  100  J  =  2,NPC 

FB  =  (  F(J.2)  +  F(J-1,2)  )/2.D+00 

UB  =  (  U(J.2)  +  U(J-1 .2)  )/2.D+00 

VB  =  (  V(J,2)  +  Vd-U)  )/2.D+00 

ZB  =  ( Z(J,2)  +  Z(J-1,2)  )/2.D+00 

PB  =  (  P(J,2)  +  P(J-1.2)  )/2.D+00 

AMB  =  (  AM(J,2)  +  AM(J-1.2)  )/2.D+00 

FVB  =  ( F(J,2)*V(J,2)  +  F(J-1.2)*V(J-1.2)  )/2.D+00 

USB  =  (  U(J,2)*U(J,2)  +  U(J-1.2)*U(J-1.2)  )/2.D+00 

FPB  =  (  F(J,2)*P(J,2)  +  F(J-1,2)*P(J-1,2)  )/2.D+00 

UZB  =  ( U(J.2)*Z(J,2)  +  U(J-1,2)*Z(J-1.2)  )/2.D+00 

DBV  =  (  B(J,2)*V(ja)  -  B(J-1,2)*V(J-1,2)  )/DETA(J  -  1) 

DBZP  =  (  BZ(J,2)*P(J,2)  -  BZ(J-1.2)*P(J-1,2)  )/DETA(J  -  1) 

IF{  NX  .EQ.  1)  THEN 

C - 

C  PREVIOUS  STATION 

C - - - 


on  nnn  nnno 


CFB  =  O.D+OO 
CUB  =  O.D+OO 
CVB  =  O.D+OO 
C:ZB  =  O.D+OO 
CPB  =  O.D+00 
CAMB  =  O.D+OO 
CFVB  =  O.D+OO 
CUSB  =  O.D+00 
CFPB  =  0.D+O0 
CUZB=  O.D+OO 
CDBV  =  O.D+00 
CDBZP  =  O.D+OO 
CRB  =  O.D+OO 
CRZB  =  O.D+OO 
ELSE 

CFB  =  ( F(J,1)  +  F(J-1.1)  )/2.D+00 

CUB  =  ( U(J.l)  +  U(J-l.l)  )/2.D+00 

CVB  =  (  V(J.l)  +  V(J-1.1)  )/2.D+00 

CZB  =  ( Z(J.l)  +  Z(J-1,1)  )/2.D+00 

CPB  =  ( P(J,1)  +  P(J-1,1)  )/2.D+00 

CAMB  =  ( AM(J,1)  +  AM(J-1.1)  )/2.D+00 

CFVB  =  ( F(J.l)*V(J.l)  +  F(J-l,l)*V(M.l)  )/2.D+00 

CUSB  =  ( U(J.1)*U(J,1)  +  U(J-1.1)*U(J-1.1)  )/2.D+00 

CFPB  =  ( F(J,l)*P(J.l)  +  F(J-1.1)*P(J-1.1)  )/2.D+00 

CUZB  =  (  U(J,l)*Z(J.l)  +  U(J-1.1)*Z(J-1,1)  )/2.D+00 

CDBV  =  ( B(J,l)»V(J.l)  -  B(J-l,l)*V(M.l)  )/DETA(M) 

CDBZP  =  (  BZ(J,l)*P(J.l)  -  BZ(J-1.1)*P(J-1.1)  )/DETA(M) 

CLB  =  CDBV  +  A1*CFVB  -  A2*CUSB  +  CAMB 
CRB  =  -CLB  +  CEL*(  CFVB  -  CUSB  ) 

CLZB  =  CDBZP  +  A1Z*CFPB 
CRZB  =  -CLZB  +  CEL*(  CFPB  -  CUZB  ) 

ENDIF 

COEFFICIENTS  OF  THE  DIFFERENCED  MOMENTUM  EQUATION  SINCE  S7(J) 
AND  S8(J)  ARE  EQUAL  TO  0.5.  THEY  ARE  SPECIFIED  WITHOUT  VECTOR  FORM 

S  i(J)  =  B(J.2)/DETA(J-1)  +  (  P1P*F(J.2)  -  CEL*CFB  )/2.D+00 
S2(J)  =  -B(J-1.2)/DETA(J-1)  +  ( P1P*F(J-1.2)  -  CEL*CFB  )/2.D+00 
S3(J)  =  (  P1P*V(J.2)  +  CEL*CVB  )/2.D+00 
S4(J)  =  (  P1P*V(J-1.2)  +  CEL*CVB  )/2.D+00 
S5(J)  =  -P2P*U(J,2) 

S6(J)  =  -P2P»U(J-1,2) 


COEFFICIENTS  OF  DIFFERENCED  ENERGY  EQUATION 

B1(J)  =  BZ(J.2)/DETA(J-1)  +  (  P1PZ*F(J.2)  -  CEL*CFB  )/2.D+00 

B2(J)  =  -BZ(J-1.2)/DETA(J-1)  +  ( PIPZ*F(J-1.2)  -  CEL*CFB  )/2.D+00 

B3(J)  =  ( P1PZ*P(J.2)  +  CEL*CPB  )/2.D+00 

B4(J)  =  (  P1PZ*P(J-1.2)  +  CEL*CPB  )/2.D+00 

B5(J)  =  CEL*(  CZB  -  Z(J.2\  )/2.D+00 

B6(J)  =  CEL*(  CZB  -  Z(J-1.2)  )/2.D+00 

B7(J)  =  -CEL*(  U(J,2)  +  CUB  )/2.D+00 

B8(J)  =  -CEL*(  U(J-1,2)  +  CUB  )/2.D+00 


DEFINITION  OF  RJ 


non 
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c - 

R(U)  =  F(J-1.2)  -  F(J,2)  +  DETA(J-1)»UB 

R(2J)  =  CRB  -  (  DBV  +  P1P*FVB  -  P2P*USB  +  CEL*(  FB*CVB-  CFB*VB  )  +  AMB  ) 

R(3J)  =  CRZB  -  (  DBZP  +  P1PZ*FPB  -CEL*(UZB+CUB*ZB-UB*CZB+CFB*PB-FB*CPB  ) ) 
R(4J-1)  =  U(J-U)  -  U(J^)  +  DETA(J-1)*VB 
R(5J-1)  =  Z(J-1.2)  -  Z02)  +  DETA(J-1)*PB 
100  CONTINUE 

57  =  5.D-01 

58  =  5.D-01 


SET  UP  DERIVED  BOUNDARY  CONDITIONS 


R(l.l)  =  0.D+00 
R(2.1)  =  O.D+00 
R(3.1)  =  O.D+00 
R(4  =  O.D+00 

R(5  J^^PC)  =  O.D+00 
RETURN 
END 

SUBROUTINE  SOLV5 

IMPLICIT  REAL*8(A-H.O-Z)aNTEGER*4(I-N) 

COMMON  /INPTl/  A1.A2A1Z^FA0^FAUT,IFLAG,GRJUA, 

1  HPLATEJlADrUS4lKVISI.RE.UINFJU 

COMMON  /ETA3/  ETAEJ4P6J«JP.DETA(200).ETA(200)A(200)J^,NNP>IXT, 

1  NX.X(600),CVW.VW(2).VWT 

COMMON  /REST4/  F(200.2).U(200.2).V(200.2)^(200;2)  J>(200,2), 

1  AM(200.2).B(200,2).BZ(200,2),BT(200) 

COMMON  /COEF5/  S  1(200), S2(200).S3(200).S4(200),S5(200).S6(200). 

1  S7.S8.B  1(200)32(200).B3(200).B4(200), 

2  B5(200)36(200).B7(200).B8(200)3(5,200) 

COMMON  /SOLV6/  DELF(200),DELU(200)  J)ELV(200)  J)ELZ(200)JDELP(200)  ABSERR 
DIMENSION  A12(200).A13(200)A14{200),A15(200), 

1  A2 1 (200),A22(2CO),A23(200).A24(200).A25(200), 

2  A3 1  (200)  A32(200).A33(200)  A34(200).  A35(200). 

3  G 1 1  (200), G 1 2(200),G  1 3(200).G  I4(200),G  1 5(200), 

4  G2 1  (200),G22(200),G23(200),G24(200),G25(200), 

5  G31(200),G32(200),G33(200),G34(200),G35(200), 

6  W1(200),W2(200).W3(200),W4(200),W5(200) 
C===================================== 

C  ELEMENTS  OF  TRIANGLE  MATRIX-ALFA  FOR  J  =  0,  SET  UP  SO  THAT 
C  DELF(1),  DELU(l),  ALFA0*DF,LZ(1)  +  ALFAl*DELO(l)  ARE  ZERO 
C  SINCE  A1 1(J)  IS  UNITY,  IT  IS  OMITTED 
A12(l)  =  O.D+00 
A13(I)  =  0.D+00 
A14(l)  =  O.D+00 
A15(I)  =  0.D+OO 
A21(I)  =  O.D+00 
A22(l)  =  1  D+00 
A23(l)  =  O.D+OO 
A24(l)  =  O.D+00 
A25(l)  =  O.D+00 
A3 1(1)  =  0.0+00 
A32(l)  =  O.D+00 
A33(l)  =  0.D+00 
A34(l)  =  ALFAO 


uuu  uuu  uuu  uuu 
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A35(l)  =  ALFAl 


ELEMENTS  OF  W- VECTOR  FOR  J  =  0 


W1(1)  =  R(1.1) 
W2(l)  =  R(2.1) 
W3(l)  =  R(3.1) 
W4(l)  =  R(4.1) 
W5(1)  =  R(5,1) 


FORWARD  SWEEP 

DO  100J  =  2J^ 

AAl  =  A(J)*A24(J-1)  -  A25(J-1) 

AA2  =  A(J)*A34(J-1)  -  A35(J-1) 

AA3  =  A(J)*A12(J-1)  -  A13(J-1) 

AA4  =  A(J)*A22(J-1)  -  A23(J-1) 

AA5  =  A(J)*A32(J-1)  -  A33(J-1) 

AA6  =  A(J)*A14(M)  -  A15(J-1) 

AA7  =  A(J)*S6(J)  -  S2(J) 

AA8  =  A(J)*S8 
AA9  =  A(J)*B6(J) 

AAIO  =  A(J)*B8(J)  -  B2(J) 

DET  =  AA4*AA2  -  AA1*AA5  -  A21(J-1)*(  AA3*  AA2  - 
1  AA5*AA6 )  +  A31(J-1)*(  AA3*AA1  -  AA4*AA6 ) 


ELEMENTS  OF  TRIANGLE  MATRIX-GAMMA 

GI  I(J)  =  (  AA5*AA1  -  AA4*AA2  +  A(J)*A(J)*(  A21(J-1)*AA2-  A31(J-1)*AA1  )  )/DET 
G12(jr)  =  (  AA3*AA2  -  AA5*AA6  -  A(J)*A(J)*(  AA2-  A31(J-1)*AA6 )  )/DET 
GISO  =  (  AA4*AA6  -  AA3*AA1  +  A(J)*A{J)*(  AAl-  A21(J-1)*AA6 )  )/DET 
G14(jr)  =  G11(J)*AI2(J-1)  +  G12(J)*A22(J-1)  +  G13(J)*A32(J-1)+  A(J) 

G15(J)  =  G11(J)*A14(J-1)  +  G12(;)*A24(J-1)  +  G13(J)*A34(J-1) 

G21(J)  =  (  S4(J)*(  AA2*AA4  -  AA1*AA5  )  +  A31(J-1)*(  AA1*AA7 
1  -  AA4*AA8  )  +  A21(J-1)*(  AA5*AA8  -  AA7*AA2 )  )/DET 

G22(J)  =  (  AA2*AA7  -  AA5*AA8  +  A31(J-ir(  AA3*AA8 
1  -  AA6*AA7  )  +  S4(J)*(  AA5*AA6  -  AA2*AA3  )  )/DET 

G23(J)  =  (  AA4*AA8  -  AAI*AA7  +  S4(J)*I  AA3*AA1 
1  -  AA4*AA6  )  +  A2I(J-1)*(  AA7*AA6  -  AA3*AA8  )  )/DET 

G24(J)  =  G21(J)*AI2(J-I)  +  G22(J)*A22(J-1)  +  G23(J)*A32(J-1)-  S6(J) 

G25(J)  =  G21(J)*AI4(J-1)  +  G22(J)*A24(J-1)  4  G23(J)*A34(J-1)-  S8 
G3I(J)  =  (  B4(J)*(  AA4*AA2  -  AA5*AA1 )  -  AA9*(  A21(J-1)*AA2 
I  -  A31(J-I)*AA1  )4- AA10*(A21(J-1)*AA5- A31(J-1)*AA4))/DET 

G32(J)  =  (  B4(J)*(  AA5*AA6  -  AA3*AA2 )  4-  AA9*(  AA2 
1  -  A3 1  (J- 1  )* AA6 )  -  AA 1 0*(  A A5-  A3 1  (J- 1  )*A^\3  )  )/DET 

G33(J)  =  (  B4(J)*(  AA3*AAl  -  AA4*AA6 )  -  AA9*(  AAl 
1  -  A21(J-1)*AA6 )  4-  AA10*(  AA4-  A2HJ-1)*AA3  )  )/DET 

G34(J)  =  G31(J)*A12(J-1)  4-  G32(J)*A22(J-1)  4-  G33(J)*A32(J-1)-  B6(J) 

G35(J)  =  G31(J)*A14(J-1)  4-  G32(J)*A24(J-1)  4-  G33(J)*A34(J-1)-  B8(J) 


ELEMENTS  OF  TRIANGLE  MATRIX-ALFA 


A12(J)  =  -A(J)-G14(J) 
A13(J)  =  A(J)*G14(J) 
A14(J)  =  -GI5(J) 


ono  ooos  ono 

II  It  O  I  • 
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A15(J)  =  A(J)*G15(J) 

A21(J)  =  S3(J) 

A22(J)  =  S5(J)  -  G24(J) 

A23(J)  =  A(J)*G24(J)  +  S1(J) 

A24(J)  =  S7  -  G25(J) 

A25(J)  =  A(J)*G25(J) 

A31(J)  =  B3(J) 

A32(J)  =  B5(J)  -  G34(J) 

A33(J)  =  A(J)*G34(J) 

A34(J)  =  B7(J)-G35(J) 

A35(J)  =  A(J)*G35(J)  +  B1(J) 

ELEMENTS  OF  W- VECTOR 

W1(J)  =  R(U)  -  G11(J)*W1(J-1)  -  G12(J)*W2(M) 

1  -  G13(J)*W3(J-1)  -  G14(J)»W4(J-1)  -  G15(J)*W5(J-1) 

W2(J)  =  R(2J)  -  G21(J)*W1(J-1)  -  G22(J)*W2(J-1) 

1  -  G23(J)*W3(J-1)  -  G24(J)*W4(J-1)  -  G25(J)*W5(J-1) 

W3(J)  =  R(3J)  -  G31(J)*W1(J-1)  -  G32(J)*W2(J-1) 

1  -  G33(J)*W3(J-1)  -  G34(J)*W4(J-1)  -  G35(J)*W5(J-1) 

W4{J)  =  R(4  J) 

W5(J)  =  R(5J) 

CONTINUE 


BACKWARD  SWEEP 

D1  =  A3l(NPC)*(  A23(NPC)*A15(NPC)  -  A13(NPC)*A25(NPC) )  + 

1  A33(NPC)*(  A25(NPQ  -  A21(NPQ*A15(NPC) )  - 

2  A35(NPC)*(  A23(NPC)  -  A21(NPC)*A13(NPC) ) 

DF  =  A35(NPC)*(  A13(NPC)*W2(NPC)  -  A23(NPC)*W1(NPC) )  +  A33(NPC)* 

1  ( A25(NPC)*  W 1  (NPC)  -  A 1 5(NPC)* W2(NPC) )  -  W3(NPC)* 

2  (  A13(NPC)*A25(NPQ  -  A23(NPC)*A15(NPC) ) 

DV  =  A31(NPC)*(  A15(NPC)*W2(NPC)  -  A25(NPC)*W1(NPC) )  -  A35(NPC)* 

1  (W2(NPC)- A21(NPC)*W1(NPQ)  +  W3(NPC)* 

2  (A25(NPC)- A15(NPC)*A21(NPC)) 

DP  =  A31(NPC)*(  A23(NPC)*W1(NPC)  -  A13(NPC)*W2(NPC) )  +  A32(NPC)* 

1  (  W2(NPQ  -  A21(NPC)*W1(NPC)  )  -  W3(NPQ*(A23(NPC)  -  A13(NPC)*A21(NPC)) 


ELEMENTS  OF  DELTA- VECTOR  FOR  J  =  NP 


DELF(NPC)  =  DF/Dl 
DELU(NPC)  =  O.D+00 
DELV(NPC)  =  DV/Dl 
DELZ(NPQ  =  O.D+00 
DELP(NPC)  =  DP/Dl 
DO200J  =  NPC-l.l.-l 

BBl  =  DELU(J-»-l)  -A(J+l)*DELV(i+l)  -  W4(J) 
BB2  =  DELZ(J-i-l)  -A(J+1)*DELP(J+1)  -  W5(J) 
CCl  =  W1(J)  -  A12(J)*BB1  -  A14(J)*BB2 
CC2  =  W2(J)  -  A22(J)*BB1  -  A24(J)*BB2 
CC3  =  W3(J)  -  A32(J)*BB1  -  A34(J)*BB2 
DDl  =  A13(J)  -  A(J+1)*A12(J) 

DD2  =  A23(J)  -  A(J+1)*A22(J) 

DD3  =  A33(J)  -  A(J+1)*A32(J) 
EEl=A15(J)-A(J-t’>*A14{J) 


65 


EE2  =  A25(J)  -  A(J+1)*A24(J) 

EE3  =  A35(J)  -  A(J+1)*A34(J) 

DETT=  DD2*EE3  +  A21(J)*DD3*EE1  +  A31(J)*DD1*EE2 
1  -  A31(J)*DD2*EE1  -  A21(J)*DD1*EE3  -  DD3*EE2 

C - 

C  ELEMENTS  OF  DELTA- VECTOR 

DELF(J)  =  (  CC1*DD2*EE3  +  CC2*DD3*EE1  +  CC3*DD1*EE2  - 
1  CC3*DD2*EE1  -  CC2*DD1*EE3  -  CC1*DD3*EE2  )/DETT 

DELV(;)  =  (  CC2*EE3  +  A2I(J)*CC3*EEl  +  A31(J)*CC1*F^2 
1  -  A3 1(J)*CC2*EE1  -  A21(J)*CC1*EE3  -  CC3*EE2  )/DETT 

DELU(J)  =  BBl  -  A(J+1)*DELV(J) 

DELP(J)  =  (  CC3*DD2  +  A21(J)*CC1*DD3  +  A31(J)*CC2*DD1 
1  -  A31(J)*CC1*DD2  -  A21(J)*CC3*DD1  -  CC2*DD3  )/DETr 

DELZ(J)  =  BB2  -  A(J+1)*DELP(J) 

200  CONTINUE 

C - 

C  NEW  VALUES  OF  F.U.V^J*  WITH  UNDER-RELAXATION 

ABSERR  =  O.D+00 

DO300J  =  i;^c 

F(J,2)  =  F(J^)  +  DELF(J)/1.2D+00 

U(J,2)  =  U(J^)  +  DELU(J)/1.2D+00 

V(J.2)  =  V(j;2)  +  DELV(J)/I.2D+00 

Z(J^)  =  Z(J^)  +  DELZ(J)/I.2D+00 

P(J.2)  =  P(J^)  +  DELP(J)/1.2D+00 

IF(  Z(J^)  .LT.  O.D+00 )  THEN 

Z(Ja)  =  O.D+00 

ENDIF 

ABSERR  =  ABSERR  +  DABS(  DEIJF(J)  )  +  DABS(  DELU(J) )  + 

1  DABS(  DELV(J) )  +  DABS(  DELZ(J) )  +  DABS(  DELP(J) ) 

300  CONTINUE 

C - - 

C  RESET  BOUNDARY  CONDITIONS  TO  AVOID  TRUNCATION  ERROR 

U(U)  =  0.D+00 

Z(l,2)  =  ( l.D+00  -  ALFA1*P(I,2)  )/ALFA0 

RETURN 

END 

SUBROUTINE  INUPD 

IMPLICIT  REAL*8(A-H.O-Z).INTEGER*4(I-N) 

COMMON  /INPTl/  Al.A2wMZ.ALFA0^FAI.IT.IFLAG,GRJUA. 

I  HPLATEJIADIUS,RKVISI,RE,UINFJ?I 

COMMON  /EOS2/  ZMF(101)3HO(101).TEP(10I)  J?HOM(200).TEMP(200), 

1  DVC(101)J’R(I01).DVCM(200).PRM(200)JC 

COMMON  /ETA3/  ETAE  W  J«’.DETA(200).ETA(200),A(200);^,NNPJ4XT, 
1  NX,X(600),CVW.VW(2).VWT 

COMMON  /REST4/  F(200.2).U(200.2).V(200,2)^(200;2)  J>(200^). 

1  AM(200.2)3(200,2).BZ(200.2).BT(200) 

C====== - === - ^====== 

IF(  IFLAG  .EQ.  0 )  THEN 

NPC  =  NP6 

ELSE 

NPC  =  NP 

ENDIF 


noon 
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SET  UP  CURVATURE  EFFECTS 


CURVAT  =  2.D+00**1.5D+00*HPLATE/RADIUS*(  X(NX)/GR  )**2.5D-01 
CURVAT  =  2.D+00**1.5D+(X)*HPLATE/RADIUS*(  XcNX)/RE  )**5.0D-01 
C  CURVAT  =  O.D+00 

RHOM(l)  =  RHO(lOl) 

TEMP(1)  =  TEP(101) 

AI.l(U)  =  l.EMX) 

BT(1)=  l.D+00 

B(1.2)  =  BT(l)*RHOM(l)*DVC(101)/RHO(l)/DVC(n 
BZ(U)  =  B(U)/PR(101) 

DO  100JC  =  2;vIPC 
CALL  SOURCE 

AM(JCJ2)  =  ( RHO(l)/RHOM(JC)  -  l.D+00 )/(  RHO(l)/ 

1  RHO(lOl)  -  l.D+00 )  •  2.D+00*X(NX)*RI 

C  AM(JC^)  =  0.D+00  INVOKE  FOR  FORCED  CONVECTION** 

BT(JC)  =  BT(JC-1)  +  CURVAT*RHO(1)/2.D+00*(  l.D+00/RHOM(JC) 

1  +  l.D+00/RHOM(JC-l)  )*(  ETA(JC)  -  ETA(JC-1) ) 

B(JC;2)  =  BT(JC)*RHOM(JC)*DVCM(JQ/RHO(I)/DVC(l) 

BZ(JC.2)  =  B(JC.2)/PRM(JC) 

100  CONTINUE 

RETURN 
END 

SUBROUTINE  SOURCE 

IMPLICIT  REAL*8(A-H,0-Z),INTEGER*4(I-N) 

COMMON  /EOS2/ZMF(101)4iHO(10l),TEP(101)JiHOM(200),TEMP(200), 
1  DVC(101),PR(101).DVCM(200)J>RM(200)JC 

COMMON  /REST4/  F(200.2).U(200.2).V(200,2)Z(200,2)  J>(200,2). 

1  AM(200, 2)3(200, 2), BZ(200,2),BT(200) 


DO  101=  1,100 

IF(  (  Z(JC,2)  .GE.  ZMF(D )  .AND.  (  Z(JC,2)  .LT.  ZMF(I+1) ) )  THEN 
RHOM(JC)  =  RHO(I)  +  (  RHO(I+l)  -  RHO(I)  )*(  Z(JC,2)  • 

1  ZMF(I))/(ZMF(I+1)-ZMF(I)) 

TEMP(JC)  =  TEP(I)  +  ( TEPa+1)  -  TEP(I)  )*(  Z(JC,2)  - 
1  ZMF(I) )/( ZMF(I+ 1 )  -  ZMF(I)  ) 

DVCM(JC)  =  DVC(I)  +  (  DVC(I+1)  -  DVC(I)  )*(  Z(JC,2)  - 
1  ZMF(I) )/(  ZMF(I+ 1 )  -  ZMF(I) ) 

PRM(JC)  =  PR(I)  +  (  PR(I+1)  -  PR(I)  )*(  Z(JC,2)  - 
1  ZMF(0  )/( ZMF(I+ 1 )  -  ZMF(I) ) 

ENDIF 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  GROWTH 

IMPLICIT  REAL*8(A-H,0-Z),INTEGER*4(I-N) 

COMMON  /ETA3/  ETAE,NP6  J4P,DETA(200)3TA(200),A(200)>JPC.NNPJ4XT, 
1  NX,X(600),CVW,VW(2).VWT 

COMMON  /REST4/  F(200,2),U(200,2),V(200,2),Z(2002)3(200,2), 

1  AM(200,2).B(200,2),BZ(200,2),BT{200) 


NNP  =  3 
NP0  =  NP 
NPl  =  NP+  1 


non 
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NP  =  NP  +  NNP 

IF(  NP  .GT.  310  )  THEN 

NP  =  400 

ENDEF 

DO  100  L=  1,2 
DO  101  J  =  NP1,NP 
DETA(J)  =  DETA(J-1) 

A(J)  =  A(M) 

ETA(J)  =  ETA(J-1)  +  DETA(J) 

F(J,L)  =  F(NP0J-)  +  U(NP0X)*(  ETA(J)  -  ETA(J  -1) ) 

U(J  J.)  =  U(NP0,L) 

V(W  =  V(NP0,L) 

Z(J«  =  Z(NP0a-) 

P(JJL)  =  P(NP0.L) 

AM(Ji)  =  AM(NP0X) 

=  B(NP0,L) 

BZ(JX)  =  BZ(NP0« 

101  CONTINUE 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  OUTPUT 

IMPLICIT  REAL*8(A-H.O-Z)  JNTEGER*4(1-N) 

COMMON  /INPTl/  A1.A2,A1ZALFA0,ALFA1,IT,1FLAG,GRJUA, 

1  HPLATEJIAD1US.RKV1S1JIE.U1NF411 

COMMON  /EOS2/  ZMF(101)JIHO(101).TEP(101)JIHOM(200),TEMP(200), 

1  DVC(101)J>R(101).DVCM(200).PRM(200)JC 

COMMON  /ETA3/  ETAE  J'IP6,NP.DETA(200).ETA(200),A(200),NPC,NNP,NXT. 
1  NX.X(600),CVW,VW(2).VWT 

COMMON  /REST4/  F(200,2),U(200,2),V(200,2),Z(200,2)  J>(200,2), 

1  AM(200,2).B(200.2).BZ(200.2).BT(200) 

DIMENSION  Y(200) 

C  - : - -,«===::;==================== 

1DNP=  1 
NDXl = 1 


OUTPUT  CONTROL  (ADD  OR  DELETE  COMMENTS  FOR  DESIRED  OUTPUT) 


1F(  (NX  .LT.  NXT)  .AND.  (NX  .NE.  1) )  THEN 
NDX  =  NDX1 
NX1  =  NX+  1 

IPRINT  =  NXl/NDX  -  NX/NDX 
1F(  IPRINT  .EQ.  0 )  GO  TO  10 
ENDIF 


C  WRITE(6,1000)  NX,X(NX),VW(2) 

1000  FORMAT(1HO,4HNX  ^.I3,5X,3HX  =JJ10.7,5X.4HVW  =£10.4) 

RMBR  =  RHOM(l)*VW(2) 

C  WRrrE(6£000)  RE£JA£MBR 

2000  FORMAT(1H0.4HRE  =£10.4,5X.4HJA  =£10.4,5X,5HMBR  =,E10.4) 

WR1TE(6£050)  X(NX)£MBR 
2050  FORMAT(lH0£10.7,T20,F15.10) 

C  WRnE(6.3000) 

3000  FORMAT(1HO,2X,1HJ,3X,3HETA.4X,1HF.9X.2HF1,8X.2HF2.8X.1HZ,9X,2HZ1) 

C  WRITE(6,4000)  ( J,ETA(J)£(J,2).U(J,2),V(J.2)7(J.2)P(J,2)4  =  IJMP.IDNP ) 
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C  WRniE(6,4050)  X(NX)JIMBR.P(U) 

4050  FORMAT(lH0J^10.7,T20.F15.10,T40J^15.10) 

C  WRITE(6,4000)  NP£TA(NP),F(NP.2).U(NP.2).V(NP,2)Z(NP,2)4"(NP,2) 

C  WRITE(6.4100)  NP£TA(NP)J^(NP^).U(NP^),V(NP.2) 

4100  FORMATa3J^.2.3Fl0.5) 

4000  FORMAT(1H043^.2,5F10.5) 

C  UPDATE  PROFILES 

C - 

PX  =  HPLATE*X(NX) 

C  WRITE(64000)  PX 

5000  FORMATaHO.PX='£15.8) 

lODO  100‘j= 

F(J.1)  =  F(J^) 

U(J.1)  =  U(J^) 

V(J.1)  =  V(J^) 

Z(J.1)  =  Z(J;2) 

P(J,1)  =  P(J^) 

B(J.1)  =  B(J,2) 

BZ(J.l)  =  BZ(J.2) 

AM(J,1)  =  AM(J^) 

C  TRANSFORM  DIMENSIONLESS  VARIABLE  TO  PHYSICAL  LENGTH 

Y(;)  =  RADIUS*(  DSQRT(  BT(J) )  -  l.D+00  ) 

UP  =  DSQRT(  4.D+00*9.80665D+00*(  RHOM(NP)/RHOM(l)  -  l.D+00 ) 

1  *X(NX)*HPLATE  )*U(J^) 

C  IF(  (NX  .NE.  1)  .AND.  (EPRINT  J9E.  0)  .OR.  (NX  .EQ.  NXT)  )  THEN 

C  IF(J£Q.  1)THEN 

C  ETAFP  =  0.D+00 

C  ELSE 

C  ETAFP  =  ETAFP  +  ( GR/X(NX)  )**2.5D-01/HPLATE/RHOM(NP) 

C  1  /2.D+00*(  RHOM(J)  +  RHOM(J-l)  )*(  Y(J)  -  Y(J-l)  )/DSQRT(  2.D+00  ) 

C  ENDIF 

C  WRITE(6,6000)  ETAFP2(J.2) 

6000  FORMAT(1H04^15.6,5XJ^15.6) 

C  ENDIF 

C  WRITE(6,7000)  Y(J)J^(J,2),UP,TEMP(J) 

7000  FORMAT(1HO,4F15.6) 

C  WRITE(6,8000)  Y(J)Z:(J^) 

8000  FORMAT(IHO,2F15.6) 

100  CONTINUE 

NX  =  NX  +  1 

C  INITIALIZE  WALL  STREAM  FUNCTION  DUE  TO  WALL  BLOWING  EFFECTS 

C  CVW  =  -HPLATE/RKVISI/DSQRT(  8.D+00  )*RHO(101)/RHO(1)/GR**2.5D-Ol 

CVW  =  -DSQRT(2.D+00)/2.D+00*(HPLATE*RHO(101)/DVC(1))/RE**0.5D+00 
F(U)  =  F(l.l)*(  X(NX-1)/X(NX)  )**5.0D-01 
1  +  CVW*(  VW(1)  +  VW(2)  )/2.D+00*(  X(NX)  -  X(NX  -  1)  )/X(NX)**5.0D-01 

RETURN 
END 


Table  A-1.  Data  for  state  relationships  used  in  flow 
calculations:  temperature  (T),  density  (p),  viscosity  (p),  and 
Prandtl  number  (Pr)  as  functions  of  conserved  scalar  (Z). 
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z 

T(K) 

p  (kg/m3) 

p  (kg/m-s) 

Pr 

0.0000 

298.2 

5.9696E-01 

1.5573E-05 

0.7484 

0.0100 

685.1 

2.5784E-01 

3.1761E-05 

0.7444 

0.0200 

1002.1 

1.7726E-01 

4.2585E-C5 

0.7466 

0.0300 

1228.2 

1.4475E-01 

4.9860E-05 

0.7448 

0.0400 

1388.1 

1.2070E-01 

5.5290E-05 

0.7605 

0.0500 

1454.7 

1.0348E-01 

5.8402E-05 

0.7864 

0.0600 

1492.2 

9.0196E-02 

6.0575E-05 

0.8088 

0.0700 

1517.7 

7.9616E-02 

6.2141E-05 

0.8258 

0.0800 

1537.2 

7.0993E-02 

6.321  lE-05 

0.8378 

0.0900 

1553.2 

6.382  lE-02 

6.3839E-05 

0.8460 

0.1000 

1567.1 

5.774  lE-02 

6.4050E-05 

0.8513 

0.1100 

1580.5 

5.2516E-02 

6.3883E-05 

0.8549 

0.1200 

1608.6 

4.8224E-02 

6.375  lE-05 

0.8616 

0.1300 

1644.6 

4.4299E-02 

6.2848E-05 

0.8664 

0.1400 

1681.9 

4.0843E-02 

6.1308E-05 

0.8674 

0.1500 

1711.6 

3.7935E-02 

5.9710E-05 

0.8643 

0.1600 

1735.9 

3.5437E-02 

5.8213E-05 

0.8583 

0.1700 

1757.2 

3.3242E-02 

5.6813E-05 

0.8507 

0.1800 

1777.1 

3.1282E-02 

5.5498E-05 

0.8418 

0.1900 

1796.4 

2.9512E-02 

5.4258E-05 

0.8321 

0.2000 

1816.2 

2.7895E-02 

5.3092E-05 

0.8219 

0.2100 

1837.5 

2.6400E-02 

5.2002E-05 

0.8113 

0.2200 

1861.9 

2.4998E-02 

5.0992E-05 

0.8006 

0.2300 

1892.1 

2.3652E-02 

5.0092E-05 

0.7879 

0.2400 

1933.3 

2.2324E-02 

4.9367E-05 

0.7790 

0.2500 

1989.1 

2.1016E-02 

4.8916E-05 

0.7689 

0.2600 

2045.0 

1.9871E-02 

4.8576E-05 

0.7597 

0.2700 

2089.0 

1.8932E-02 

4.8117E-05 

0.7510 

0.2800 

2125.2 

1.8114E-02 

4.7582E-05 

0.7426 

0.2900 

2159.3 

1.7371E-02 

4.7034E-05 

0.7343 

0.3000 

2193.1 

1.6674E-02 

4.6523E-05 

0.7262 

0.3100 

2229.8 

1.6005E-02 

4.6095E-05 

0.7186 

0.3200 

2273.0 

1.5342E-02 

4.5806E-05 

0.7114 

0.3300 

2330.0 

1.4657E-02 

4.5769E-05 

0.7048 

0.3400 

2415.6 

1.3895E-02 

4.6242E-05 

0.6997 

0.3500 

2541.7 

1.3048E-02 

4.7426E-05 

0.6972 

0.3600 

2685.4 

1.2242E-02 

4.8919E-05 

0.6969 

0.3700 

2827.4 

1.1538E-02 

5.037  lE-05 

0.6977 

0.3800 

2961.2 

1.0939E-02 

5.1719E-05 

0.6988 

0.3900 

3062.0 

1.0499E-02 

5.2666E-05 

0.7002 

0.4000 

3111.1 

1.0236E-02 

5.2969E-05 

0.7016 

0.4100 

3119.2 

1.0095E-02 

5.275  lE-05 

0.7030 

0.4200 

3097.5 

1.0040E-02 

5.2164E-05 

0.7042 

0.4300 

3053.9 

1.0048E-02 

5.1308E-05 

0.7052 

0.4400 

2997.4 

l.OlOOE-02 

5.0313E-05 

0.7058 
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Table  A-1  --  continued 


0.4500 

2938.2 

1.0170E-02 

4.9332E-05 

0.7063 

0.4600 

2884.1 

1.0236E-02 

4.8482E-05 

0.7069 

0.4700 

2838.8 

1.0288E-02 

4.7819E-05 

0.7079 

0.4800 

2802.1 

1.0321E-02 

4.7338E-05 

0.7092 

0.4900 

2772.6 

1.0340E-02 

4.7007E-05 

0.7109 

0.5000 

2748.4 

1.0346E-02 

4.6788E-05 

0.7129 

0.5100 

2728.1 

1.0344E-02 

4.6648E-05 

0.7152 

0.5200 

2710.5 

1.0336E-02 

4.6565E-05 

0.7176 

0.5300 

2694.8 

1.0324E-02 

4.6520E-05 

0.7201 

0.5400 

2680.3 

1.0309E-02 

4.6499E-05 

0.7227 

0.5500 

2666.6 

1.0293E-02 

4.6495E-05 

0.7254 

0.5600 

2653.2 

1.0276E-02 

4.6497E-05 

0.7281 

0.5700 

2639.7 

1.0261E-02 

4.6496E-05 

0.7309 

0.5800 

2625.6 

1.0247E-02 

4.648  lE-05 

0.7336 

0.5900 

2610.5 

1.0238E-02 

4.6442E-05 

0.7362 

0.6000 

2593.5 

1.0234E-02 

4.6363E-05 

0.7387 

0.6100 

2573.4 

1.0240E-02 

4.6216E-05 

0.7411 

0.6200 

2548.0 

1.0263E-02 

4.5959E-05 

0.7433 

0.6300 

2513.0 

1.0316E-02 

4.5504E-05 

0.7451 

0.64C0 

2459.7 

1.0431E-02 

4.467  lE-05 

0.7461 

0.6500 

2375.5 

1.0658E-02 

4.3217E-05 

0.7460 

06600 

2267.5 

1.0997E-02 

4.1296E-05 

0.7452 

0.6700 

2159.3 

1.1377E-02 

3.9400E-05 

0.7448 

0.6800 

2063.1 

1.1753E-02 

3.7799E-05 

0.7454 

0.6900 

1982.8 

1.2096E-02 

3.6584E-05 

0.7474 

0.7000 

1917.8 

1.2395E-02 

3.574  lE-05 

0.7506 

0.7100 

1864.9 

1.2654E-02 

3.5193E-05 

0.7550 

0.7200 

1820.7 

1.2881E-02 

3.4857E-05 

0.7601 

0.7300 

1782.7 

1.3086E-02 

3.4669E-05 

0.7657 

0.7400 

1748.9 

1.3275E-02 

3.458  lE-05 

0.7718 

0.7500 

1718.1 

1.3454E-02 

3.4563E-05 

0.7782 

0.7600 

1689.3 

1.3628E-02 

3.4590E-05 

0.7848 

0.7700 

1662.0 

1.3800E-02 

3.4645E-05 

0.7916 

0.7800 

1635.5 

1.3972E-02 

3.4713E-05 

0.7985 

0.7900 

1609.6 

1.4148E-02 

3.478  lE-05 

0.8056 

0.8000 

1583.8 

1.4331E-02 

3.4838E-05 

0.8127 

0.8100 

1573.2 

1.4384E-02 

3.4304E-05 

0.8110 

0.8200 

1562.1 

1.4442E-02 

3.3771E-05 

0.8094 

0.8300 

1550.1 

1.451  lE-02 

3.3258E-05 

0.8083 

0.8400 

1539.0 

1.4588E-02 

3.2700E-05 

0.8061 

0.8500 

1532.9 

1.4657E-02 

3.1937E-05 

0.7995 

0.8600 

1526.4 

1.4729E-02 

3.1148E-05 

0.7927 

0.8700 

1519.6 

1.4806E-02 

3.0333E-05 

0.7856 

0.8800 

1512.3 

1.4888E-02 

2.949  lE-05 

0.7782 

0.8900 

1504.4 

1.4976E-02 

2.8623E-05 

0.7706 

0.9000 

1496.0 

1.5070E-02 

2.7728E-05 

0.7627 

0.9100 

1486.8 

1.5172E-02 

2.6807E-05 

0.''547 

0.9200 

1476.8 

1.5284E-02 

2.5862E-05 

0.7467 

0.9300 

1465.9 

1.5406E-02 

2.4893E-05 

0.7385 

0.9400 

1453.7 

1.5542E-02 

2.3906E-05 

0.7305 

0.9500 

1440.2 

1.5696E-02 

2.2904E-05 

0.7225 

0.9600 

1425.0 

1.5872E-02 

2.1896E-05 

0.7149 

Table  A- 1  --  continued 


0.9700 

1407.8 

1.6075E-02 

2.0894E-05 

0.7078 

0.9800 

1388.2 

1.6315E-02 

1.9914E-05 

0.7014 

0.99^ 

1365.8 

1.6602E-02 

1.8979E-05 

0.6960 

1.0000 

1340.5 

1.6946E-02 

1.8120E-05 

0.6920 
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APPENDIX  B 

SUFHCIENT  AND  NECESSARY  CONDITION  FOR  SIMILARITY  SOLUTIONS 

To  obtain  a  sufficient  and  necessary  condition  for  similarity  solutions,  the 
nonsimilar  terms  of  the  right-hand  sides  of  the  momentum  and  mixture  fraction 
equations,  Eqs.  (2.30)  and  (2.38),  are  examined.  First,  the  nonsimilar  term  in  Eq.  (2.30) 
is  set  equal  to  zero. 

^f-^r  =  0  (B.l) 

34 

The  dimensionless  stream  function,  f,  is  assumed  to  consist  of  two  pans  using  separation 
of  variables. 

f=X(4)H(ri)  (B.2) 

Substituting  into  Eq.  (B.l),  one  obtains 


XX'(H’2  -  HH")  =  0 


(B.3) 


Three  possibilities  exist  for  Eq.  (B.3),  X  =  0,  X’  =  0,  or  H'2  -  HH"  =  0.  The  condition  of 
X  =  0  leads  to  trivial  solution  and  is  excluded  from  funher  consideration.  The  solution  to 
H'2  -  HH"  =  0  is  obtained: 
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H'(0)  =  0  H'(«)  =  1  (B.5) 

Equation  (B.4)  subject  to  the  boundary  conditions  of  Eq.  (B.5)  also  leads  to  a  trivial 
solution. 

Finally,  X'  =  0  is  satisfied  if  X  =  constant  or  Vw  is  proportional  to  xhe 

dimensionless  stream  function  is  a  function  of  T)  only,  implying  similarity  solutions  for 
momentum  equation. 

To  complete  the  necessary  condition,  the  nonsimilar  term  of  Eq.  (2.38)  is  set  to 
zero.  Separation  of  variables  for  mixture  fraction  is  employed. 

—  f-  — Z'  =  0  (B.6) 


Z  =  Xi(^)Hi(ti) 


(B.7) 


